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ABSTRACT 

Factorization is an approach to linear programming (LP) in which the algebraic 
elements of the LP tableau are organized in such a way that a large portion of the 
tableau may be represented implicitly and generated from the remaining explicit 
part. In dynamic row factorization, the row structure of the LP model instance 
influences the algebraic structure of the tableau, and the dimension of the algebraic 
elements may change as the solution progresses. 

We present three algorithms motivated by this approach, each resulting from 
a different LP model row structure: generalized upper bound (GUB) rows, pure net- 
work rows and generalized network rows. We describe implementations of all three 
algorithms, specifying data structures for tableau and basis inverse representations 
and detailing procedures for manipulation and update of these representations. 

Computational results are presented for a number of real-world models taken 
from a variety of applications and industries. From each model, one or more partic- 
ular instances are solved by each of our three implementations and by a commercial- 
quality mathematical programming system. The characteristics of the four solvers 
are compared and contrasted. Previous research on related algorithms by others 
suggests that these algorithms are properly viewed as specialized approaches, useful 
only on narrow classes of problems. Our computational results strongly refute this 
view, and instead suggest that each algorithm is superior to the general simplex 
approach on a wide range of problem classes and structures. 
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I. INTRODUCTION 



A. INTRODUCTION 

A recurring theme in the development of simplex-based algorithms for linear 
programming has been the identification and exploitation of special problem struc- 
ture. Ideas as apparently disparate as the simplex method for bounded variables, 
primal and dual decomposition methods, pure and generalized network primal sim- 
plex algorithms and primal partitioning schemes may be unified to a degree by 
interpreting their development in this context. 

The factorization approach due to Graves and McBride [1976] provides a unify- 
ing framework from which we may reinterpret many existing algorithms and, through 
the application of the common principles embodied in the approach, develop new 
algorithms. This approach has as its central thesis the idea of recognizing, isolat- 
ing and exploiting special structures which may occur in a certain type of linear 
programming tableau. An example of such special structure is generalized network 
rows. Generalized network rows are naturally specified as a row structure in which 
each column has at most two nonzero elements within the rows. 

This paper adopts the design principles and algorithmic structure suggested 
by Graves and McBride [1976] to develop algorithms that, while general in nature in 
the sense that each may be used to solve any linear programming problem instance, 
are strongly tailored to exploit a particular row structure. We call this approach 
“dynamic row factorization"; "row factorization” because we exploit the structure 
in the basic tableau which is induced by the row structure of the LP model instance, 
and “dynamic” because the dimension of the structure may vary (or even fail to be 
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present) as the solution progresses. In our setting, we require the row structure of 
the model instance to be specified prior to solving, and to remain fixed throughout 
the solution process. An extension of this approach is to allow the row structure of 
the model instance to vary as the problem is solved. While we do not develop an 
implementation of the second approach, we show that it is a conceptually simple 
extension of our work. 

Each algorithm is developed by factoring the constraints of the LP model 
into two classes: those that have special structure (factored) and those that do 
not (explicit). This factoring of constraints induces a factored structure in the LP 
tableaus which may be exploited computationally. We treat each of three structures 
of factored constraints in this work: generalized upper bounds (GUB), pure networks 
and generalized networks. 

We implement each of the three algorithms by integrating it within the struc- 
ture of a state-of-the-art mathematical programming system: the X-System of 
Brown and Graves [1975]. We do so both to demonstrate the feasibility of imple- 
menting such methods and to determine whether or not such methods have practical 
value for the practitioner interested in solving real-world problems. Commercial im- 
plementations of similar GUB algorithms have generally failed to inspire enthusiasm 
and are apparently falling into disuse [Kennington [1978]]. While the first commer- 
cial implementation of a related “pure network with side constraints” algorithm is 
(to our knowledge) just now becoming available, reports of research implementations 
have generally supported the view that such algorithms have limited application. 
The reports of implementations of similar “generalized network with side constraint” 
algorithms have offered less promise than their pure network counterparts. 

After we review the related literature, we provide a detailed presentation of the 
underlying tableau-based algorithm which forms the foundation of all our subsequent 
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work. We do so because the algorithm is not widely known and may be unfamiliar 
to the reader. We then review the general factorization approach of Graves and 
McBride [1976]. We provide a design template for our developmental approach, and 
then present the dynamic row factorization algorithms for GUB, pure network and 
generalized network row structures. Computational results are presented, and we 
then summarize our conclusions and suggest avenues for further research. 

B. SURVEY OF RELATED LITERATURE 

While the terms “partitioning” and “factorization” are frequently used inter- 
changeably in the literature, we observe a distinction between the two approaches. 
We consider partitioning methods to be those that are based on special structure 
in the original problem instance. This structure in the problem instance need not 
induce special structure into the LP tableau, and in fact the method need not be 
tableau-based. This is in contrast to our view of factorization, in which the algo- 
rithm is based on special structure which occurs in the simplex basis and thus in the 
basic tableau. We thus classify Dantzig- Wolfe [1960] and Benders [1962] Decompo- 
sitions and Rosen's [1964] primal partitioning method as examples of partitioning 
methods. 

Perhaps the first example of what we consider factorization is the treatment of 
simple upper bounds by Dantzig [1954] and, independently, by Charnes and Lemke 
[1954]. They observe that it is more efficient to enforce the “structural” simple 
upper bound constraints with logical tests within the algorithm rather than treat 
them explicitly along with other constraints. While not originally presented in the 
context of a formal tableau factorization, the approach is easily viewed as such and 
is consistent with the general approach. 
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The mutual primal-dual method of Graves [1965] focuses attention on the 
special role of nonnegativity constraints in linear programming. A clear distinction is 
drawn between the computational convenience of treating nonnegativity constraints 
implicitly rather than explicitly and the unambiguous mathematical equivalence 
of all problem constraints, structural or nonnegativity. Emphasizing the special 
importance of inequality constraints, the approach yields an elegant theory (see 
Graves [1987]) and, as we will see, efficient implementations. We view this algorithm 
as the first formal example of factorization. 

Dantzig and Van Slyke [1967] extend the factorization approach applied earlier 
to simple upper bounds in a more structured manner in their treatment of general- 
ized upper bounds (GUB). In a problem with p GUB constraints and m structural 
constraints, their approach requires a working basis inverse of dimension (m + 1), a 
considerable savings when p is large. 

Hartman and Lasdon [1972] specialized this approach to the multicommodity 
capacitated transshipment problem. In this case, the structure of the basic pure 
network columns introduces additional structure into the working basis, allowing 
further simplifications in basis representation and update techniques. Helgason and 
Kennington [1977] develop techniques for representing the working basis inverse in 
product form and provide graphic interpretation of the graph updates. Kennington 
[1977] reports an implementation of the algorithm. 

McBride [1972] and Graves and McBride [1976] formalize and generalize the 
factorization approach. They view it is as a unifying framework for tableau-based 
simplex specializations and illustrate this by developing a variation of the GUB 
algorithm of Dantzig and Van Slyke [1967] and a GUB algorithm for doubly coupled 
linear programs of Hartman and Lasdon [1970]. They present a new algorithm 
for the set partitioning linear programming problem and an equality-constrained 
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form of the pure network with side constraints model. McBride [1972] reports an 
implementation of a GUB factorization. 

Schrage [1975] extends the approach of simple and generalized upper bounds 
by considering variable upper bounds (VUB), which are constraints of the form 
Xj < Xk, where x* is said to be the variable upper bound of xj. His algorithm allows 
the implicit representation of the VUB constraints by expressing VUB variables in 
terms of other basic columns. This permits the basis representation to be treated 
in two parts, one part a large matrix which changes infrequently and thus needs 
to be updated only occasionally, and the second part a small working basis which 
requires regular attention. Thus, computation and storage savings may be realized. 
Schrage [1978] extends these ideas to what he calls generalized variable upper bounds 
(GVUB) constraints, which arise frequently in models involving fixed charges. 

Klingman and Russell [1975] sketch a factorization method for solving trans- 
portation problems with side constraints. They suggest techniques for performing 
simplex iterations and updating the problem representation. Chen and Saigal [1977] 
present a similar approach for solving capacitated network flow problems with ad- 
ditional linear constraints. Both the above presentations consider a graph-theoretic 
view of the basis update mechanism and allow the basis representation to be treated 
in two parts, a part which corresponds to a rooted spanning tree defined on the un- 
derlying graph, and a general working basis inverse. Glover, Karney, Klingman and 
Russell [1978] report an implementation of the Klingman and Russell design, but 
one which (curiously) only accommodates a single side constraint. McBride [1989] 
reports an implementation which requires the pure network rows to be equalities 
and allows more than one side constraint. 
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The problem of generalized network problems with side constraints is ad- 
dressed by Hultz and Klingman [1976]. They present details for the simplex price- 
out, column generation and basis update. Hultz and Klingman [1978] report an 
implementation that (curiously) solves the “singularly constrained” generalized net- 
work problem. McBride [1989] reports an implementation that is not restricted to 
a single side constraint. 

The factorization approach has been extended by the consideration of em- 
bedded structures. Glover and Klingman [1981] consider a linear program which 
contains embedded pure network structure, i.e., the pure network structure appears 
in only a subset of the rows and columns of the technological coefficient matrix of 
the problem. Their approach yields an algorithm similar in spirit to the algorithms 
for the pure network with side constraint model, but the presence of the “side vari- 
ables” significantly complicates the basis representation and update. They report 
an implementation of the algorithm but (curiously) restrict the problem suite to 
problems having no complicating variables. 

McBride [1985] treats the problem of linear programs with embedded general- 
ized network structure. He presents methods for pricing, column generation, basis 
representation update and data structures. A successful implementation is reported 
which is approximately five times faster than MINOS [1977] for the models tested. 

Interest in developing algorithms to solve problems with special substructures 
has been accompanied by work to identify such substructures in problem instances. 
Greenberg and Rarick [1974] and Brown and Thomen [1980] develop algorithms 
to identify GUB sets. Brown and Wright [1984] develop algorithms for identifying 
pure network constraint substructures. Brown, McBride and Wood [1985] present 
a method for locating generalized network structures, both embedded and row-only 
structures. 
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Todd [1983] examines factorization from a geometric standpoint and constructs 
a geometric interpretation which is in large measure equivalent to the algebraic 
development of Graves and McBride [1976]. 



II. MUTUAL PRIMAL-DUAL METHOD 

A. INTRODUCTION 

This chapter presents the mutual primal-dual linear programming method in- 
troduced by Graves [1965] which provides the algorithmic framework and notational 
conventions for the research which follows. 

We begin with a complete algebraic development of the primal algorithm fol- 
lowed by a less detailed symmetric discussion of the corresponding dual algorithm. 
We conclude with a unified treatment of the two which establishes the theoretical 
importance of the algorithm and justifies its use as the foundation of our specializa- 
tions. 

The following presentation scrupulously restates the Graves [1965] algorithm, 
incorporates later discussion by McBride [1972], and accommodates large-scale im- 
plementations by Brown and Graves [1975]. The view presented here is not available 
in standard reference texts, and is included in the interest of completeness. 

B. PRIMAL PROBLEM STATEMENT (PLP) 

The traditional statement of the linear programming (LP) problem is: 

(LP) min : wy 

s.t. a,y < r, , i = 1, . . . ,m 

£jV> 0 , j = l,...,n, 

where y is an n-vector of decision variables, w is an n-vector of cost coefficients, 
each a, is an n-vector of technological transformation coefficients, each r, is a scalar 
right-hand side coefficient, and e 3 is the j ih unit vector. While this statement of 
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the problem is clear and unambiguous, there are reasons for preferring an alter- 
native. The insistence upon drawing a formal distinction between the “structural” 
constraints a,y < r,- and the “nonnegativity” constraints ejy > 0 obscures the math- 
ematical structure of the problem by suggesting that the two types of constraints 
are inherently different. Certainly the exploitation of the special structure of the 
ejy > 0 constraints leads to computational efficiencies in the implementation of the 
algorithm. However, in the theoretical development of the algorithm, we prefer to 
treat them simply as general inequality constraints. 

In order to achieve a consistent form, we rewrite the nonnegativity constraints 
35 — e^y < 0 and group them with the structural constraints. The problem statement 
then becomes: 

(PLP) min : wy 

s.t. : a.iy<ri , i = 1, . . . , m + n, 

where wy is called the extremal function. The constraints of (PLP) define a set of 
feasible solutions which we shall call the feasible set, F: 

F = {?/ € R n | a { y < r,-, i = 1 , . . . , m + n } . 

Since F is the intersection of a finite number of closed half-spaces, F is itself a closed 
set. If F is nonempty and bounded, then an optimal solution to (PLP) will occur 
at an extreme point of F. 

A point y° € F is said to be a feasible point or feasible solution. If, for 
constraint i, a,y° < r, . constraint i is satisfied at y° , and the quantity r, — a,y° is 
the slack in constraint i at y° . If, on the other hand, for constraint t, a t y° > r, , 
constraint i is violated at the point y° , the magnitude of the violation being a,y° — r, 
(the negative of the slack). 
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A point y° E R r ‘ is defined to be a basic solution of (PLP) if there exists an 
independent subset {a tl , a,- 2 , . . . , a,„} of {aj, a 2 , . . . ,a m+n } such that a,v y° = r tj for 
j = 1, . . . n. Such an independent subset {a tl ,a,- 2 , . . . , a,„} of {aj, a 2 , . . . ,a m+n } is a 
basis for R n , and y° is the (basic) solution to the system a tj y° = r tj , j = 1, . . . , n. 

Each basic solution of (PLP) corresponds to an extreme point of F and for each 
extreme point of F there is at least one corresponding basic solution of (PLP). Since 
there are at most ( m n n ) ways of choosing an independent subset of n vectors from 
{ai,a 2 , .. . ,a m + n }, the number of basic solutions of (PLP) and thus the number of 
extreme points of F is finite. Hence, (PLP) can be solved by searching among its 
basic solutions. The algorithm to be developed here will implicitly enumerate the 
set of basic solutions of (PLP) and terminate in one of three states: 

1. F = 0 

(no feasible solution exists) ; 

2. the extremum is unbounded 

(for every real number a there is a point y° € F w'ith wy° < a ); or 

3. there exists at least one optimal solution 

(a point y’ 6 F with ivy’ < wy V y € F ). 

We will first consider the problem of finding a basic feasible solution to (PLP). 
Having achieved this, we will then consider the task of finding a basic feasible solution 
which is also optimal. 

C. OBTAINING FEASIBILITY 

Since we have included the nonnegativity constraints — ejy <0,j = l,...,n 
in our structural constraints a t y < r, , i = 1, . . . . m -fi n, and since the origin is the 
unique solution to the independent system: 
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-ejy = 0, j = l,...,n, 



the origin is a basic solution that is always immediately available for (PLP). 

Let y° be any basic solution to (PLP). By definition, there are at least n linearly 
independent constraints which are exactly binding at y° . Among the m remaining 
constraints, typically some will be satisfied at y° and others will be violated. Our 
strategy will be to focus on one violated constraint at a time, which we will call 
the target constraint. Moving from one basic solution to another, we will attempt 
to reduce the violation of the target constraint until it becomes satisfied. We will 
restrict our choices of basic solutions in such a way that all constraints that are 
already satisfied at y° remain satisfied at each subsequent basic solution. Once the 
target constraint becomes satisfied, we then select some other violated constraint 
as the new target constraint, and repeat the process. We proceed until either all 
constraints are satisfied and we have obtained a basic feasible solution, or we find 
a violated constraint which cannot be satisfied, in which case we conclude that no 
basic feasible solution exists. 

To formalize these ideas, define S(y°) to be the set of indices of all constraints 
satisfied at a basic solution y° : 



S{y°) = {l < f < m + n | a,y° < r,} . 

Of course, | S(y°) | > n. 

Suppose constraint k is violated at the basic solution y°. Then aky° > r* and 
k S(y°). A necessary condition for the existence of a basic feasible solution is that 
there exists a basic solution y with: 
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%) 3 S(y°) 



( 2 . 1 ) 



and 



akV < a k y°. (2.2) 

If there exists no such basic solution y satisfying (Eq. 2.1) and (Eq. 2.2 ) we conclude 
that F = 0. Thus, we may restrict our attention to basic solutions satisfying (2.1) 
and (2.2). 

Let {a,,, a,- 2 , . . . ,<z, n } be a basis for R n at y°. For notational convenience we 
will partition the constraints into two sets, those that are basic at y° and those that 
are nonbasic at y°: 



' h ■ 

b 2 


_ 


1 

' ** ft ft 

to " ' 

1 


; / = 


i 

^ ^ ... 
i 


— 


’ r n ' 


.K. 




. <*.•„ . 




1 

. < 

i 




1 

e 

i 



D = 


’ <h 1 
d 2 


= 


1 

fl *n+ 2 


; 9 = 


i — ' — 

... 




' r.n+l 
r <„+2 


(2.4) 




£ 

^3 




CU , 

L l n|m J 




. Qm _ 




v ■ 

*n+m - 





Define B and T> to be the set of indices of the rows of B (basic constraints) 
and D (nonbasic constraints) respectively. Using this notation, the current basic 
solution y° may be expressed as By 0 = /, and since the rows of B are by definition 
linearly independent, B~ l exists and y° = B~ x f . 

Let y be any other point of R n . Then y can be expressed as: 



y = y° + {y - y 0 )- 
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We have chosen this representation since at y° the basic constraints are satisfied 
exactly and thus it is the direction vector ( y — y°) of any proposed move from y° 
to y that determines whether or not a given basic constraint will remain satisfied 
at y. ( Of course, it is the magnitude of (y — y°) that is important in determining 
whether the nonbasic constraints remain satisfied as well.) 

Now, for a given arbitrary basis {p’,p 2 ,. . . ,p"} of R n and any point y, there 
must exist scalars Aj, . . . , A n such that: 



y = y° + (y - y°) = y° + H = y° + P\ . 

Since the choice of P is arbitrary, let us choose P = — B -1 , which we call the 
conjugate row basis. Then: 



By = B(y° + P A) = B(y° - B’ 1 A) = By 0 - A = / - A, 

and thus y satisfies the basic constraints By < f if and only if A > 0. 

Thus, the set of points in R n which satisfy the basic constraints 6,-y < /; for 
i € B can be characterized as: 



yeR n \y = y o + Y, x :F’ X :^ 0 ’j£ B 



with P = -iT 1 . 

It then follows that every point y satisfying 



S(y) 2 5(y°), 



(2.5) 
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can be expressed as 



y = y° + PA, 



for some A > 0 and P = — B~ l . 

To address condition (2.2), assume that we have selected constraint k as our 
target constraint (and thus k is violated at y° , i.e. , d k y° > g k ). To reduce the 
violation of the target constraint, we seek a point y such that: 

d k y < d k y°, (2.6) 

holds. Since such a point must also satisfy (2.1), it must have a representation of 
the form y = y° + PA for some A > 0. Replacing y by y° + PA , we have: 

d k (y° + PA) < d k y°, 

which holds if and only if d k P A < 0. Since A > 0, it follows that a necessary 
condition for the existence of a point y satisfying (2.5) and (2.6) is that at least 
one element of the vector d k P be negative, i.e., d k pd < 0 for some j, 1 < j < n. If 
d k P > 0, we conclude that F = 0. 

Suppose that d k p l < 0. Then (2.6) holds at all points of the form: 

y = y° + a ip 1 , A I > 0. 

The generic point y° + A ip 1 lies along an edge of the set of all points satisfying the 
constraints which are basic at y° . If there is more than one / for which d k p l < 0, 
any one can be chosen. For the purpose of exposition, we will designate the first 
such index encountered as l 7". 
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The violation in constraint k decreases linearly as Aj increases, 
becomes satisfied when: 



Constraint k 



dk(y° + A ip ') = g k , 

or when A* = t , with 



t A (g k - d k y°) 

dkP 1 

A geometric interpretation is that y = y° -j- tp 1 is the point at which the ray 
jy € R n | y = y° + A /p ( , A; > oj pierces the hyperplane {y € R n \ d k y = g *}. 

Define A* to be our ultimate choice for A / . Choosing t as the value for Af 
satisfies the target constraint k, but we are required by (2.5) to continue to satisfy 
all nonbasic constraints which are already satisfied at y° (if any). The choice of t 
may cause violations of such constraints. Thus, A' must also satisfy the condition: 



(2.7) 



<f,(y° + V)<& v jfDnsi,"). 



Writing this as 



VjV < 9 : ~ d : y°, 



we find that we must choose A / < s, where 



A 

s = min 
T>nS(y° 



liiLlMl | djp' > o} . (2.8) 

'S(y °) \ d 3 p l J 

If this set is empty, define s — oc . Of course, it is possible for s to be equal to zero. 
Thus, by choosing 
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A; = min {t, 5 } , (2.9) 

where t and s are as defined in (2.7) and (2.8) respectively, we obtain the largest 
reduction in the violation of the target constraint consistent with the feasibility 
restriction (2.5) with respect to the chosen direction p l . 

The selection A; according to (2.9) leads to a new basic solution y 1 = y 0 +A*p , r 
If A; > 0 at each step, the transitivity of (2.5) and (2.6) guarantees that no basic 
solution will be repeated before any given target constraint is satisfied, or until 
the conclusion is reached that F = 0. At any given step, either the feasible set is 
shown to be empty, a basic solution is found which satisfies the target constraint 
or a positive reduction in the violation of the target constraint is achieved as the 
result of moving to another basic solution. Since the set of basic solutions is finite, 
the third alternative may be repeated a finite number of times before one of the 
first two alternatives occurs. If this constraint becomes satisfied, (2.5) guarantees 
that every subsequent solution will also satisfy the constraint. We may then select 
a new violated nonbasic constraint as the target constraint. Since there are a finite 
number of constraints, we will either discover a basic feasible solution or determine 
in a finite number of steps that none exists. 

If A* = 0 at any step of the algorithm, we say that a block has occurred. 
Blocks will be discussed later in detail. 

D. FROM FEASIBILITY TO OPTIMALITY 

Once feasibility is achieved, say at y°, the process of proceeding to optimality 
can be thought of as minimizing over F the violation of the constraint xvy < wy° — a , 
where a is a sufficiently large positive constant. 
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Let y° be a basic feasible solution with By 0 = / and P = —B~ l . Since every 
point y £ R n which satisfies By < f may be written as y = y° + P\ for some A > 0, 
the value of the extremal function wy at such a point y is: 



wy = w(y° + P A) = wy° -f wPX = wy° + ^ A jWp 3 . 

j = i 

Thus, a necessary condition for the existence of a point y £ F such that wy < wy° 
is: 



wp 3 < 0 for some j, 1 < j < n. 

If wp 3 > 0, then y° must be a feasible, optimal solution to (PLP). 

Suppose wp 1 < 0. Since all constraints are satisfied at y°, the greatest reduc- 
tion in the value of the extremal function in the direction p 1 is achieved when A* is 
chosen to be s, where s is as defined in (2.8). If $ = oo , (PLP) has an unbounded 
extremum. If 0 < A* < oo , a positive reduction in the value of the extremal function 
can be achieved by moving from the basic feasible solution y° to the basic feasible 
solution y = y° -f X'p 1 . If A* = 0, a block is encountered. 

Once a feasible solution has been obtained. (2.5) and (2.6) become equivalent 
to: 



»EF , (2.10) 

wy < wy° . (2.11) 

The transitivity of (2.10) and (2.1 1) ensures that no basic feasible solution will 
be repeated as long as A” > 0 at each step. At each iteration, either a basic feasible 
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solution is found to be optimal, the solution is found to be unbounded or a positive 
reduction in the value of the extremal function is achieved by moving to another 
basic feasible solution. Since the set of basic feasible solutions is finite and no basic 
feasible solution is repeated, the third alternative can only occur a finite number of 
times. Thus, if A* > 0 at each step and if a finite optimal solution exists, it will be 
found in a finite number of steps. 

E. BLOCKS 

Suppose that y° is a basic solution satisfying 



By 0 = f, 

where, as before, B is of dimension n by n and nonsingular and p 7 is the j th column 
of P = — B~ l . Then, y° is the unique point in R n lying in the intersection of the n 
hyperplanes: 



{y e R n \ b;y° = fi, i = l,...,n}. 

Suppose that p 7 is a direction that either leads to a reduction in the violation 
of some target constraint or, if y° is feasible, a reduction in the value of the extremal 
function. Since the nonbasic constraints are of the form: 



d,y° <9,-, i € V, 



a block occurs when, for at least one nonbasic constraint k, 



dkV° = 9 k- 



IS 



and 



d k p‘ > 0. (2.12) 

Thus, any movement away from y° in the direction p l will result in violation of the 
nonbasic constraint d k y° < g k . 

If d k y° = g k , then y° is in a sense “over-determined”, and y° is said to be 
a “degenerate” solution. Geometrically, y° lies in the intersection of at least n + 1 
hyperplanes. Algebraically, there is more than one basis that can be formed from 
the row vectors { ai,fl2 , • • • , a m +n} for which: 



a i, V° = r i} , j = l,...,n, 

holds, and there is more than one basic solution which corresponds to the extreme 
point y° . A block is therefore encountered when an over-determined solution sat- 
isfies (2.12) in an “improving” direction (e.g., one that leads either to a reduction 
in infeasibilitv or. if feasible, to a reduction in the value of the extremal function). 
We will develop a method for dealing with blocks later. 

F. THE PRIMAL ALGORITHM AND A SUPPORTING BASIC TAB- 
LEAU 

The primal algorithm proceeds as follows: 

1. Identify an initial basic solution. Notice that the origin satisfies 

-/• y° = o, 



and thus may always serve as the initial solution. 
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2. If the current solution is infeasible, select a violated primal constraint index k 
(i.e., an index k for which g k - d k y° < 0). This requires the quantity: 

g-Dy° . 

3. Within the target row, select an element of the proper sign (negative, according 
to our convention) whose index, /, specifies a “transformation column”. If the 
current solution is infeasible, this requires the quantity: 

d k P , 

with the transformation column satisfying d k p l < 0. If the current solution is 
feasible, this requires the quantity: 



wP 

with the transformation column satisfying wp l < 0. If no such element exists, 
then the problem is infeasible (if the current solution is infeasible and d k P > 0) 
or optimal (if the current solution is feasible and wP > 0). The task of 
selecting such an index / is commonly referred to as “pricing” or a “pricing 
strategy”. 

4. Compute t as in (Equation 2.7). If the current solution is feasible, assign 
t = oc. 

5. Compute s as in (Equation 2.8). This computation is commonly called a “ratio 
test' , or a “minimum ratio test”. This computation requires the quantities: 

g — Dy° and Dp 1 . 

G. Compute as in (Equation 2.9). If \j = oc, the problem is unbounded. 
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7. Update the current solution according to the computation: 

y 1 = y° + A 

8. Update the assignments of constraint indices to B and T>, and update P. 

9. Go to Step 2. 

The only quantities needed at each step of the algorithm are the matrix DP , 
the column vector g — Dy° and the row vector wP. For notational convenience, we 
can define D to be an (m + 1) by n matrix whose bottom row d m+ i is w, and g to 
be a (m + 1) dimension column vector whose bottom entry, ^ m+1 , is zero. Then 



* i o c, o o 

gm +1 ~ d m+ iy = 0 -wy = -wy , 

which is the negative of the extremal function value at the point y° . We can now 
conveniently display all the relevant information by forming the (m + 1) by (m + 1) 
matrix: 



DP | g-Dy], 



(2.13) 



which we will call the basic tableau. 

By computing the basic tableau in partitioned matrix form, we may isolate the 
important algebraic components required by this method. To do so. let us assume 
that at the current basic solution the basis consists of h structural constraints and 
(n — h) nonnegativity constraints. Then (2.3) can be written as: 
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and / = 



' b , 






b 2 




di 2 


h 




a 'h 


bh+ 1 






bh + 2 













‘/i 




r u ' 


A 






A 




r .h 


fh+1 




0 


fh+2 




0 


Jn 




_ 0 



In partitioned matrix form, 



and thus 



B = 



An 

0 




A 12 
-I 



\ 

/ 



}h 

} n — h , 



n — h 






P = -B~ l = 



—A~ l 
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— An An 



\ o 

Similarly, (2.4) can be written as: 



}h 

} n — h 





' di 




t J1 




9 1 




' 0 




A 




-e J2 




9 2 




0 


D = 


d h 

j 




~ e jh 


and g = 






0 




dh+ 1 








5/l+l 




r < h+ l 




dh+? 




a 'h + 2 




5/1+2 




r *h + 2 




I 

ft- * 

3 

i 




d l n-f m 




Qm 




r »n + m 



and in partitioned matrix form: 



(2.14) 
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D = 



Then 



- 1 

\ ^21 



\ 

0 

^22 J 



}h 

} m — h 



( 2 . 15 ) 



DP = —DB~ X = 






, 4-1 
'mi 



n — h 

A u x A 12 

- 1 A A /t-1 



\ 



}h 

} m — h , 



\ — -^ 21-^11 -^22 — A 2 \A xi A \2 } 

and we shall call DP the principal part of the tableau. By partitioning to = 
(101,102)) g T = (<7i)<?2) T and y° = (yi,y%), basic tableau ( 2 . 13 ) may be written 
in partitioned matrix form as: 



n—h 






A -1 

'mi 



y? 



All A12 

— A21AU A22 — A21A11 A12 g2 — -^21^1 — A22V2 



\ 



-W1A11 1 w 2 - W1A11A12 



-wy 



J 



}h 

} m — h 

}l ■ 



( 2 . 16 ) 



Note that 1/° is displayed explicitly in the tableau ( 2 . 16 ). Also, t/° = 0 since the 
corresponding nonnegativity constraints are basic and thus binding. 

As the algorithm proceeds from one basic solution to another, ( 2 . 16 ) can be 
updated using a slight variant of the Gauss-Jordan transformation (which we will 
call a pivot). The transformation is as follows: 



Let 



DP\g — Dy° = [o,j] be the basic tableau at any basic solution y° . If 



v s i = d s p l 0 then the basic tableau at the basic solution: 



1 0 1 ( ^s.n+1 \ l 0 , 

y = y + \ -^-y- Jp = y + 



d.v 1 ) 
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written as 




is calculated as: 



general elements ( i ^ s,j ^ l) : 

, v 3J v,j 

”■> = - — ■ 

pivotal column (i ^ s, j = l) : 



-v,i 

v u = — > 

v >,i 



pivotal row ( i = s, j ^ /) : 



v s j 



pivotal element i = s, j = l 

1 



v s,l = 



V a ,l 



After applying the transformation to (2.16), row and column interchanges may 
be necessary to restore the explicit structure of B in (2.14) and D in (2.15). A proof 
of this may be found in Graves [1987], who also extends the basic algorithm to 
include free variables, equality constraints and bounded variables. 



G. DUAL LINEAR PROBLEM (DLP) 

Corresponding to (PLP) is the following linear program called its dual: 

(DLP) max : xr 

X 

s.t. : xA < w 

xl < 0. 

The relationship between (DLP) and (PLP) is as follows: 
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1. if (PLP) is infeasible, then (DLP) is either infeasible or unbounded. 



2. if (PLP) is unbounded, then (DLP) is infeasible. 

3. if (PLP) possesses a finite optimal solution y‘ , then (DLP) also possesses a 
finite optimal solution x* and wy' = x'r. 

See Graves [1987] for a proof of this. 

Solving (DLP) directly has important advantages in certain problems, and we 
will show that its solution plays an important role in dealing with blocks in (PLP). 
We will thus develop a dual algorithm in a manner similar to that for (PLP), sparing 
some of the symmetrical detail where possible. 

Let x 0 be a basic solution for (DLP). Then there exists an independent subset 



{a J1 , a 12 , . . . , a? m } of {a 1 , a 2 , . . . , a m+n } such that x 0 a?' = w ^' , i = l,...,m, and 
this linearly independent subset will be called the basis for R m at x 0 . 



Partitioning the dual constraints into those which are basic at x a and those 
which are nonbasic at x 0 yields: 




( 2 . 17 ) 



and 




( 2 . 18 ) 



Define T and K to be the set of indices of the columns of T (basic dual 
constraints) and I\ (nonbasic dual constraints), respectively. The current basic 
solution x 0 may then be expressed as 
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x Q = xlT 1 , 



and any other point x in R m may be expressed as 



X = Xo + (x — Xo ). 

For a given arbitrary basis {91,92? • • • , 9 m} of R m , any such point x may be repre- 
sented as: 



m 

X = Xo + ^ 9 ,-^' = Xo + %1>Q. 

t=i 

The current basic constraints are x 0 a J * = tiF* , i = and for these con- 

straints to remain satisfied at the generic point x = x Q -f ibQ, we must have 
[xa J ’’ < vji ' , i = 1, . . . , 7tj]. Choosing Q = T -1 as a convenient basis for R m at x 0 , we 
then have: 



xT = (x 0 + tpQ)T = x 0 T + tpQT = u + ij>, 



and thus we must have ib < 0. 

Now consider a violated dual nonbasic constraint t , where a violation means 



x 0 k l > v l . 

To reduce the violation, we seek a point x at which 



xk\ < x 0 k J . 



2G 



holds. Such a point must also satisfy the basic constraints, and thus the condition 
for a reduction in the violation of the dual target constraint is: 



xk l = ( x 0 + ipQ)k l = x 0 k l + ip Qk l < x 0 k\ 



which implies that 



ipQk 1 < 0. 

Since ip < 0, a necessary condition for constraint t to be satisfied is that there exists 
an i with > 0. If none of the elements of q{k 1 are positive, we conclude that the 
dual constraints are inconsistent. 

The dual algorithm proceeds along the edges 



x = x 0 + ip'q ,, 

from basic solution to basic solution. We define 



S{x a ) = {l<j<m + n| x 0 a j < ur 3 } , 

and insist that 



S(x) D S(x 0 ). 



To satisfy this condition, we consider the effect of moving along the edge of a 
general nonbasic dual constraint A' 3 which is currently satisfied at x 0 . We must have 
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xk 3 = (x 0 + x l>'qi)k : = x 0 k 3 + xp'q { k 3 < v 3 . 
Since the constraint k 3 is satisfied at x 0 , we have 



x 0 k 3 < v 3 =$■ v 3 — x a k 3 > 0, 

and since t/>* < 0, if qik 3 < 0 then as xj>' decreases in value we approach the boundary 
of the constraint. If v 3 - x 0 k 3 = 0 then q { k 3 < 0 causes an immediate violation if 
xj)' < 0. Thus, a block has occurred. 

We have shown that our choice for denoted xpl, should be 

V’* = max {6, c} , (2.19) 

where 

b = {v ( - Xotfyqik*, (2.20) 

and 

c = max {(t^ — x 0 k 3 )/qik 3 \ v 3 — x 0 k 3 > 0, q;k 3 < o) . (2.21) 

Once dual feasibility has been achieved, we wish to maximize: 



xr = x a r + xp(Qr). 

Since V’ < 0, a necessary condition to achieve an increase in the value of the extremal 
function is that there exists an i such that q,r < 0. Hence, when <?,r > 0, we conclude 
that x a is dual optimal. 

The dual algorithm then proceeds as follows: 
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1. Identify an initial basic solution. Notice that the origin satisfies 

x° • I = 0, 

and thus may always serve as the initial solution. 

2. If the current solution is infeasible, select a violated dual constraint index t 
(i.e., an index t for which v l — Xok* < 0). This requires the quantity: 

v — x 0 K. 

Column t is then referred to as the “target column”. If the current solution is 
feasible, the right-hand side column is designated the “target column”. 

3. Within the target column, select an element of the proper sign (positive, ac- 
cording to our convention) whose index, f, specifies a “transformation row”. 
If the current solution is infeasible, this requires the quantity: 

Qk* , 

with the transformation row satisfying g,A* ( > 0. If the current solution is 
feasible, this requires the quantity: 

Qr , 

with the transformation row satisfying <?,r > 0. If no such element exists, then 
the problem is infeasible (if the current solution is infeasible and q{k l < 0) or 
optimal (if the current solution is feasible and q,r < 0). The task of selecting 
such an index i is commonly referred to as “pricing”, or a “pricing strategy”. 

4. Compute k as in (Equation 2.20). If the current solution is feasible, assign 
b = — oc. 
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5. Compute c as in (Equation 2.21). This computation is commonly called a 
“ratio test”. This computation requires the quantities: 

v — x 0 K and Qk 1 . 

6. Compute rpl as in (Equation 2.19). If ipl = oo, the problem is unbounded. 

7. Update the current solution according to the computation: 

*1 = zo + rp\qi. 

8. Update the assignments of constraint indices to T and K . , and update Q. 

9. Go to Step 2. 

This development shows that at each step, the quantities QK , ( v — x a K) and 
Qr are needed to proceed with the dual algorithm. 

To develop a mat rix partitioned form of the dual tableau, we proceed as before. 
Assuming the dual basis consists of h structural constraints and m — h nonnegativity 
constraints, we have: 



so 



T =(<>, /*,...,/"*) 
u = (u 1 , u 2 , . . . , u m ) 



= (a JI , a -’ 2 , . . . , a* h , e ' h+1 , e , ' ,+2 , . . . , e ,m ) 

= (w n ,w 3h ,0,0, — ,0), 



u = (it 1 ,u 2 ) = (te^O). 

The nonbasic constraints are then: 

AT = (A - 1 . k\ . . . , k n ) = (e 11 , e’ 2 , . . . , e* h , a* h+i , a- 7h+2 , . . . , a Jn ) 
r = (r’.r 2 ?-•") = (0, 0, .... 0, , tt^ +J ..... u Jn ), 
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SO V = (v 1 , v 2 ) = ( 0 ,u> 2 ) . 

The matrix partitioned form of the basis is then 



h m — h 



T = 



mi 



0 



A 2 i I ] 

and with the choice of Q = T -1 



}h , 

} m — h , 



h m — h 




\-A 2 iA^ / / }m — h , 



with the remaining constraints forming 



I< = 



I A12 
0 A 22 



The principal part of the dual tableau is: 



QI< 



A 



-1 

11 



— AnA^ 



A 



-1 

11 





■ / 


1 

CN 




1 

0 


1 

CN 

CN 



A-y' A X 2 



— A 2 \Ayy A 22 ~ A 2 \Ayy A\2 






which we find to be exactly the principal part of the primal tableau, so 



DP = QK . 



The quantities (i’ — x 0 I\) are 



v — x 0 I\ = v — (uQ)K 



31 



The quantity Qr is: 



— / 111 1 ~Au A\2 

0 / 



v-uDP 
(0,u; 2 ) 
t>P — uZ)P 
(u-uD)P 
j(0,u> 2 ) - (to'.O) 
{(0,«, 2 )-(-f',0)}P 
wP. 



— uDP 



-I 0 
A 2 I A 22 



Qr = g 



r l 

r 2 



= Q< 

= Q 4 # + 



r i 

0 



+ 



0 

r 2 



/ v4i2 
0 j422 

= Q {$ + #/} 



= Qg + QKf 



An 0 

— A2lA 1 i 1 



0 

r 2 



+ DPf 



= g + D(Pf) 
-9- D y 0 - 



0 

r 2 



0 

r 2 



+ PP/ 



We thus find that the quantities required for the dual tableau (QK, v—x 0 K and 
Qr) are exactly the quantities required for the primal tableau (DP, wP and g — Dy°). 
Therefore, the primal and dual algorithms use the same tableau. Reviewing the 
summaries of the primal and dual algorithms, we see the strong symmetry between 
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the two. Operating on a single tableau, the primal algorithm identifies a target 
row, selects a transformation column, performs a ratio test, and performs a pivot 
update, while the dual algorithm identifies a target column, selects a transformation 
row, performs a ratio test and performs a pivot update. This remarkable symmetry 
allows us not only to perform either algorithm on a single tableau, but to switch 
between one algorithm and the other as often as we choose. This fact has important 
implications for our implementation and it enables us to deal with blocking in a 
systematic and consistent way. 

H. RESOLUTION OF BLOCKING 

The resolution of blocking in either the primal or dual algorithm will be accom- 
plished by shifting to the alternate algorithm when blocking occurs. The alternate 
algorithm is applied to a subproblem of the original and at worst we are lead to a 
contracting sequence of problems to which we alternately apply the primal and dual 
algorithm. A strict contraction is assured, and thus in at most a finite number of 
steps resolution is achieved. 

We will demonstrate this procedure by assuming that we have started with 
the primal algorithm. When a block occurs, assume we have rearranged the rows 
so that the zeros in the right-hand side column occur contiguously. 

Let Aq = n + 1 and A 2 = t w'here t is the index of the target row (if the current 
solution is primal feasible, A 2 = m + 1, the extremal function row). Let A 3 be the 
pivot column which caused the block. The situation is shown in Figure 2.1. 

Define Subproblem (1) as consisting of all the columns of the original problem, 
but only those rows with zeros in the right-hand side column, and as shown in Figure 
2.2. “extremal function 1 ’ row A 2 , the target row of the original primal problem. Row 
A '4 is the row containing the blocking element. Now apply the dual algorithm to 
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Figure 2.1: Primal Block in Column k 3 



k 3 k t 

k 4 






Subproblem (1). Because the right-hand side column of any potential pivot row is 
zero, the right-hand side column of the original primal problem is invariant to any 
pivot in Subproblem (1). Let v x] denote the element in row i and column j of the 
tableau. We proceed with the dual algorithm on Subproblem (1) until one of the 
following situations occur (which may be after zero or more dual pivots): 

1. Ujt 2 ,* 3 > 0. A pivot in column k 3 no longer reduces the violation of primal 
constraint k 2 , and thus the primal block has been eliminated. We thus return 
to the original primal problem and seek a new column in which to pivot. 

2. Vk<,k 3 < 0. A pivot in column k 3 no longer threatens to violate primal con- 
straint k 4 . We thus compute new ratios using columns k j and k 3 of the original 
primal tableau to determine if we can make a gain on the target constraint k 2 
using column k 3 . 
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Figure 2.2: Subproblem (1) 
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3. We encounter a block in the dual algorithm applied to Subproblem (1). This 
occurs when a column of Subproblem (1) contains a negative element in row 
k 4 and a zero in row k 2 , as illustrated in Figure 2.3. 

Let k s be the column in Subproblem (1) causing the dual block. We now 
define a further contraction, Subproblem (2), and switch to the primal algorithm. 
Subproblem (2) consists of all the rows of Subproblem (1), but only those columns 
of Subproblem (1) with zeros in its target row (k 2 ). Note that the current target 
column in Subproblem (1) (k 3 ) must have a negative element in its “extremaF row 
(^2). As a notational convenience we reverse the sign of this element and record this 
action by labeling the column as -k 3 . Column -k 3 becomes the right-hand side 
column of Subproblem (2), as shown in Figure 2.4. 

The number of columns in Subproblem (2) is reduced by at least one from 
the number of columns in Subproblem (1). We now apply the primal algorithm to 
Subproblem (2) until (which may be after zero or more primal pivots): 
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Figure 2.3: Dual Block in Subproblem (1) 
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Figure 2.4: Subproblem (2) 
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1. v* 4i _* 3 > 0, implying that u/t 4i * 3 < 0. A pivot in column k 3 no longer threatens 
to violate primal constraint k 4 . We have thus removed the original primal 
block, and we proceed by returning to the original primal problem, and com- 
puting new ratios using columns ki and k$ . 

2. vjt 4 ,A 5 > 0. A pivot in row k 4 no longer threatens to violate dual constraint 
k$, and thus we have removed the dual block in Subproblem (1). We revert to 
Subproblem (1) and continue with the dual algorithm. 

3. We encounter a block in the primal algorithm applied to Subproblem (2). We 
illustrate this situation in Figure 2.5. 

Let ke be the row causing the primal block in Subproblem (2). Since the the 
target row in Subproblem (2) (k 2 ) is identically zero for all but the right-hand side 
column ( — £3), primal unboundedness cannot occur. Also, the entire k 2 row in the 
original problem is invariant to pivots in this subproblem. We proceed just as we did 
when faced with this situation in the original primal problem. The current target 
row ( k 4 ) in the current problem (Subproblem (2)) becomes the extremal function 
row in a newly defined contraction (Subproblem (3)). All columns of the current 
problem are retained, but only those rows with zeros in the right-hand side column 
of the (— k 2 ) current problem are carried forward to the contraction. This again 
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Figure 2.5: Primal Block in Subproblem (2) 
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Figure 2.6: Subproblem 3 



assures a strict reduction in the number of rows. The new subproblem is shown in 
Figure 2.6. We then proceed with the dual algorithm, exactly as before. 

The construction of the contracting subproblems through the nested sets of 
zeros in the columns and rows guarantees a monotonic decrease in the sizes of the 
higher-order subproblems. This ensures the ultimate resolution of degeneracy and 
gives us a complete, symmetric and finite algorithm for the solution of LP problems. 

The blocking resolution scheme given here is a constructive algorithm to iden- 
tify strictly improved solutions. The restricted subproblems ultimately yield a pivot 
sequence satisfying all higher-order criteria. Geometrically, we systematically search 
the degenerate subspace for an improved representation. This is in sharp contrast 
to ad hoc “anti-degeneracy” and “anti-cycling” schemes which invoke arbitrary sec- 
ondary mechanisms not at all related to the geometry or mathematics of the problem 
at hand, and consequently admit nuisance degenerate pivots with no constructive 
motivation. 



I. RELATIONSHIP BETWEEN PRIMAL-DUAL ALGORITHM AND 
SIMPLEX METHOD 

We now depart from the presentation of Graves [1965] to discuss the relation- 
ship between the Mutual Primal-Dual Method and the classical Simplex Method. 
After a brief discussion of the similarities, we will explain our reasons for adapting 
the Primal-Dual view. 
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The Simplex Method assumes primal feasibility, and we must identify a start- 
ing basic feasible solution. Because this is seldom practical, the usual approach is to 
formulate and solve a new linear program closely related to the original which has 
the same optimal solution (assuming one exists) and possesses an easily identified 
basic feasible solution. This related problem is derived from (PLP) by augmenting 
it with slack variables, resulting in the following: 



(APLP) 


min : wy + 
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Is = r 


} m 
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In the classical Simplex view, a basis is a collection of m linearly independent 
columns. Let Ab be such a simplex basis which corresponds to a primal feasi- 
ble solution for (APLP). Suppose we partition the coefficient matrix in a manner 
determined by the basic slack variables, yielding (after possible row and column 
interchanges): 



where 



[A \ I] = 



An Aj2 

An A22 



h 0 

0 I 2 , 



Ab = 



An 

A21 



0 

h , 



and the matrix of nonbasic columns is 



Aa'b 



/] An 
0 A 22 , 



3S 



The basic variables are j \b — (yi,s 2 ) while the nonbasic variables are ysB — 
(si,y 2 ), and the principal part of the Simplex tableau is \A-b • 

Borrowing from the perspective of the Primal-Dual algorithm, we may generate 
the Simplex tableau by partitions: 
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and the tableau becomes 
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which is precisely the principal part of the Primal-Dual tableau. This is no surprise 
when we note that the basis for the Dual Algorithm, T, is exactly Ab , and the Primal 
and Dual Algorithms share the same tableau. We may interpret the Primal and Dual 
Algorithms as simply two different perspectives of the same tableau, wherein the 
Primal Algorithm a pivot is viewed as exchanging primal constraints and in the 
Dual Algorithm a pivot is viewed as exchanging dual constraints. The classical 
Simplex Method may then be interpreted as solving (PLP) using the 
Dual Algorithm perspective. That the classical Simplex Method is naturally 
interpreted as a dual algorithm comes as a surprise to the conventionally trained. 
However, the consequent mathematical insight is compelling, especially in light of 
the notational simplification and apparent underlying role of Af, 1 . 

There are several reasons for preferring the Primal- Dual Algorithm to the 
Simplex Method. From a computational standpoint, because slack variables are 
carried logically rather than introduced explicitly, we are able to clearly identify the 
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essential information needed to execute the algorithm. The matrix A plays a key 
role in the calculation of the tableau, and the entire tableau can be constructed from 
Aii and original problem data. Since Aj'/ is a submatrix of Ag 1 , it is smaller and 
requires fewer arithmetic operations to update than does Ag 1 . 

A second advantage of the Primal-Dual Algorithm lies in the flexibility it 
offers for specialization to particular problem classes or structures. Indeed, it is the 
special structure and simplicity of the nonnegativity constraints that motivate the 
development of the algorithm in the first place. It is frequently the case that other 
special structures can be identified in classes of (PLP). Examples of such structures 
include simple upper bounds, generalized upper bounds, variable upper bounds, 
pure and generalized network substructures, etc. Such structure may be static in 
that its nature and dimension remains fixed throughout the solution process, or 
the structure may be dynamic in which case its precise nature and/or dimension 
may vary as the problem is solved. Some special structures may be more strongly 
characterized by their column structure and others by their row structure. The 
Primal-Dual Algorithm allows one to effectively exploit virtually any such problem 
structure in a natural manner and greatly simplifies the implementation of such a 
specialization. 

When the linear programming problem appears as a subproblem in a more 
sophisticated solution setting (for example, in a mixed integer programming problem 
or a nonlinear programming problem), the row/column symmetry of the Primal- 
Dual Algorithm is of critical importance in specializing the solution approach. The 
inherent symmetry of the algorithm permits easy adaptation to branch-and-bound 
and cutting-plane approaches to mixed integer programming, to column generation 
settings, as well as to primal and dual decomposition techniques. 
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We believe the reason for this flexibility offered by the algorithm lies in its 
more complete mathematical foundation. There is a natural consistency that 
arises from the choice of a vector space having the same dimension as the 
problem variables that is lacking in other approaches. A natural geometric 
interpretation of the solution trajectory follows directly from this development. In- 
cidental issues such as finding an initial basic feasible solution and dealing 
with degeneracy are resolved constructively in this mathematical frame- 
work. Other approaches resort to unnecessarily complicated tangential efforts. 

All the research results reported here can be developed, with some effort, in 
the framework of the classical Simplex Method. However, we choose to present these 
results in the manner of their development - the mutual Primal-Dual view presented 
by Graves. 
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III. IMPLEMENTATION DESIGN OVERVIEW 



A. INTRODUCTION 

We seek to demonstrate the efficiency of row factorization (in particular, using 
dynamic-dimension bases). Accordingly, we will implement the ideas developed here 
and extensively test them within a commercial-quality optimization system: the X- 
System [1975] employs the Graves mutual primal-dual algorithm in a variety of 
large scale optimization applications, including linear, nonlinear, mixed integer and 
decomposed models. The benchmark test suite is drawn from a wide variety of actual 
applications, and our goal is to improve the efficiency of an already well-known and 
highly regarded system. 

B. DESIGN CONSIDERATIONS 

We want to test our ideas by repeatedly solving many medium- to large-size lin- 
ear programming problems (i.e., approximately 8,000-10,000 constraints and 15,000- 
20,000 structural variables). Larger problems are of interest, but for purposes of this 
research, we are limited to a relatively modest computer budget on an IBM 3033/AP 
computer under the VM/CMS operating system using VS Fortran 1.4.1. We wish 
to support the computational enhancements common in commercial mathematical 
programming systems (e.g., bounded variables, ranged constraints, parametric pro- 
gramming, etc.). We require a primal-dual implementation that offers complete 
flexibility in determining solution strategy. In addition, each experimental imple- 
mentation must support all the routine housekeeping of the optimization system 
(e.g., re-inversion, crash starts, relaxations, restrictions, extrinsic and intrinsic enu- 
meration. etc.). Finally, we seek the capability to identify desired row factorization 
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structure within the LP model instance, either by communication from the modeler 
or automatically. 



C. DESIGN TEMPLATE 

To establish a conceptual framework for the evaluation of our algorithms, it 
is useful to outline the important aspects of our implementation and identify the 
crucial steps which most strongly influence performance. Recall from Chapter 2 our 
basic tableau in partitioned form: 
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where we have made the substitution y° = B 1 f in the right-hand side column: 



y° = B~ l f =► (y°,y 2 ) T = 

=>(yly°2) T = (ATfo, of. 



A\\ A11 A 12 



-1 



r 1 

0 



The data requirements of the algorithm are: 

1. Access to the original problem data; 

2. A representation of that part of the current tableau we have chosen to represent 
implicitly, and: 

3. A representation of that part of the current tableau we have chosen to represent 
explicitly. 



We now consider each of these requirements in greater detail. 
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Primal simplex implementations typically are column-oriented and thus re- 
quire column-wise access to the problem data. In the primal-dual method, restric- 
tion to either column-wise or row-wise access alone exacts a serious computational 
penalty. Thus, our design allows both column- wise and row- wise access. 

The design of an efficient large-scale tableau element generation representation 
is remarkably complex. Many subtle software engineering and hardware environment 
issues can have a profound influence on the intricacy and performance of promising 
designs. The design excursions reported here are inexorably influenced by architec- 
ture of the host computer, operating system and implementation language. However, 
the reported design philosophy has been tempered with experience on many other 
computers of widely varied designs. The proposed innovations adapt quite well to 
floating-point pipelines, large cache memories and parallel architectures. 

Because of its fundamental role in the construction of the tableau, we maintain 
an explicit representation of A n (we thus refer to yin as the “explicit kernel”). 
Although any of the various techniques such as LU, LDL T or QR decomposition 
or product form inverse (e.g., Golub and Van Loan [1985] or Magnanti [1976]), are 
suitable for representing A n or Afj 1 , a difficulty arises in this algorithmic setting. 
While in the primal Simplex method all updates to the basis take the form of rank 
one column exchanges, our setting admits more general updates, which include single 
row exchanges, single row and column exchanges, single row and column deletions, 
and multiple row exchanges of A ^ (multiple column exchanges of An). While 
this does not preclude the use of any particular representation, it adds a level of 
complexity not usually encountered in more traditional implementations. 

We also maintain an explicit representation of the right-hand side column and 
the bottom (cost) row. Because of the symmetric nature of the mutual primal-dual 
method, a sensible approach is to allocate a single storage array for both quantities, 
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which taken together are called the “problem rim”. If we generate an explicit rep- 
resentation of the pivot row and pivot column of the tableau at each iteration, then 
we may update the problem rim array using the simple pivot transformation. By 
adopting the convention of labeling the nonbasic constraints as rows 1 through m 
and labeling the basic constraints as rows (m + 1) through (m + n) (assuming our 
problem instance has m structural constraints and n nonnegativity constraints), the 
problem rim array is partitioned as follows: 

1. The first portion of the array holds values corresponding to the nonbasic con- 
straints in region (?) of the tableau (3.1). Since these are nonnegativity con- 
straints and they are nonbinding (nonbasic), the values in region (?) of the rim 
array are those of the currently (possibly) nonzero variables. 

2. The second portion of the array holds values corresponding to the nonbasic 
structural constraints in region (?’?'), and thus the values in the rim array are 
the current slack or violation in these constraints. 

3. The third region of the rim array, beginning at position m -f 1, corresponds 
to basic (binding) structural constraints, and thus the values are those of the 
corresponding (possibly) nonzero dual variables. 

4. The fourth region of the rim array corresponds to basic nonnegativity con- 
straints, and thus the values are those of the corresponding (possibly) nonzero 
dual variables, conventionally called “reduced costs”. 

The rest of the tableau we represent implicitly, by simply recording in a stor- 
age array the current ordering of nonbasic and basic constraints. When additional 



information from the tableau (for example, a row or column) is required, we con- 
struct it from our representation of /4J', 1 , the current row ordering and the original 
problem data. 

An overview of the solution process is as follows: 

1. Identify an initial basic solution. As stated previously, the origin is always 
such a solution, and thus we may always begin with: 



if = -1 
If = A, 

or any other suitable basis. 

2. Check for optimality. If there arc currently no primal violations ( a negative 
value in the firs! or second region of the rim array) and there arc no dual 
violations (a positive value in the third or fourth region of the rim array), 
then the current solution is optimal. Otherwise, we proceed to step 3. 

3. Select either a primal or a dual violation, perform a pivot which makes progress 
towards reducing that violation and return to step 2. 

Since we desire to maintain current information in the rim array by means of 
the pivot transformation updates, we require a representation of the pivot row 
and pivot column at each iteration. Since we explicitly maintain only and 
the problem rim, we see that a key computational step in our implementation 
is the generation of tableau rows and columns. 

Recall from Chapter 2 the principal part of the primal tableau (a symmetric 
development from the dual perspective is also possible, and is of course equivalent): 
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where P = — B~ l is the conjugate row basis, 
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and D contains the nonbasic rows 
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D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 

Now consider the generation of column s from (3.2). Rewriting (3.2) in a 
manner that highlights our intentions: 
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( 3 . 5 ) 



^ — ^2i ( i 1 ) — An(Aii An) + Ail ) • 

By properly sequencing our computations we will exploit the fact that region 

(??) of a given column is simply a linear combination of region (?) of the same column. 

Assume we want to place the current representation of column s into a work 

array r. which we partition as : T = (r],c 2 ) r to correspond to (3.5). Expressed in 

terms of an explicit transformation kernel .4 7 /. we first compute region (?) of column 



s as: 



= MJiT 



if s is in (j). 
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or 



Zi = A^(A u ) a if 5 is in (jj). 

Having done this, we then compute z 2 as: 



Z 2 = -v4 2 l2i 



if s is in (j), 



or 



z 2 = -A21Z1 + (A 22 ) s if 5 is in (jj). 

Assuming an LU representation of An, we first compute region (t) of column 



s as: 



or 



LU Zi = e 3 if column s is in (j), 



LUzi = (A ]2 ) 5 if column s is in (jj), 

Having done this, we then compute z 2 as: 



or 



z 2 — — A 2 \Z\ 



if column s is in (j), 



z 2 = —A 2 \Z\ -fi (A 22 )’ if column s is in (jj). 

Then the current representation of column s is available in z T = (2i,2 2 ) r . 
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E. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 

The computation of row t of the tableau proceeds in a similar manner. We 
now view the principal part of the tableau as: 
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If we want to place the current representation of row t in a work array z 
partitioned as z = (23,24), we may first compute region ( j ) of row t as: 



^3 — (Ai 1 )* 



if row t is in (t), 



or 



z 3 = (—A 2 i)tA 11 1 if row t is in (u). 

We then compute: 



^4 — z 3 ( A \ 2 ) 



if row t is in (i), 



or 



Z4 = 23(^12) + (A2 2 ) t if row t is in (n). 

Alternately using an LU representation of A u , 
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z 3 LU — e t 



if row t is in (t), 



or 



z 3 LU = (— v4 2 i)t if row t is in (ii). 

We then compute z 4 as: 



z 4 = 23(^12) if row t is in (i), 



or 



z 4 = 23(^12) + (A 22 )* if row t is in (it), 

and the current representation of row t is available in z — (23, z 4 ). 

We see that in each case calculations proceed by first using a representation of 
A n to compute a portion of the row or column and then using this initial computa- 
tion and original problem data to compute the remaining part. We will discover that 
our specializations extend this approach by introducing additional tableau partitions 
which allow this computational strategy to be applied on a larger scale. 

As previously mentioned, an important implementation challenge in this algo- 
rithmic setting is the dynamic behavior of A n . We see from the primal row basis 
(3.3) and nonbasic rows (3.4) that the dimension of A\\ corresponds to the number 
of basic structural constraints, or, equivalently, to the number of nonbasic nonneg- 
ativity constraints (recall that if a nonnegativity constraint is nonbasic and thus 
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nonbinding, the corresponding variable may possibly be nonzero). Recalling that 
our primal view of a pivot is as an exchange of constraints between B and D , we 
see that one of four cases may occur during a pivot: 

1. A structural constraint enters the basis B and a structural constraint leaves 
the basis and enters D. Since the number of basic structural constraints (and 
the number of nonbasic nonnegativity constraints) remains unchanged, the 
dimension of An is unchanged. A pivot of this type involves a row in region 
(j) of B (3.3) and a row in region (it) of D (3.4), and thus it corresponds to 
a pivot coordinate in the location ((ii), ( j )) of the tableau (3.5). 

2. A nonnegativity constraint enters the basis and a nonnegativity constraint 
leaves the basis. Again, the dimension of A u remains unchanged. Since this 
pivot involves a row in region (jj) of (3.3) and a row in region (i) of (3.4), the 
corresponding tableau (3.4) pivot coordinate lies in ((?'), (jj)). 

3. A structural constraint enters the basis and a nonnegativity constraint leaves 
the basis, and thus the number of basic structural constraints (equivalently, 
the number of nonbasic nonnegativity constraints) increases by one. The di- 
mension of A n is increased by one. This corresponds to a pivot coordinate in 
region ((H), (jj)) of the tableau (3.5). 

4. A nonnegativity constraint enters the basis and a structural constraint leaves 
the basis, and thus the dimension of An is decreased by one. The correspond- 
ing pivot coordinate in (3.5) is ((?),(;))• 

We see that we may exert some influence on the behavior of the dimension of 
A u by our strategy for selecting target violations for primal and dual constraints 
(i.e.. our pricing strategy) and through our tie-breaking rules for choosing pivot 
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row/columns, and that this dynamism is an inherent feature of our algorithm. We 



have already seen the fundamental importance of the kernel (An) in our compu- 
tations. Thus, a successful implementation must manage this dynamic behavior 
efficiently and reliably. 

To illustrate in a familiar setting the challenge this offers, consider a LU rep- 
resentation of the kernel An' 



A n = LU. 



(3.6) 



A pivot in tableau coordinate ((?'), ( jj )) results in a column exchange in A\\. Writing 
(3.6) in the more convenient form: 



L~ l A n = U. 



(3.7) 



When column a k replaces the p ih column of A u to form A u , we have 



\ 



* u 




\ 



where a = L~ l a k has replaced the p ih column of U. We now must restore the 
upper triangular form of U. We would prefer a method which demonstrates strong 
numerical stability and which also preserves the sparsity of U. Several methods 



have been proposed (Bartels and Golub [1969], Forrest and Tomlin [1972], Saunders 
[1976] and Reid [1982]). Perhaps the most widely used method is that of Saunders. 

Assume for the moment that A\\ has a block (or “bump”) triangular form (we 
discuss how this is done shortly). Then A u appears as shown in Figure 3.1. 




Figure 3.1: Bump Triangular Form of A n 

Each bump consists of (possibly) several columns that extend above the main 
diagonal. These columns are called “spikes", and the tallest spike within a bump is 
placed in the right-most column, so that a bump appears as shown in Figure 3.2. 

The block triangular form of A\\ results in a LU decomposition which has the 
form shown in Figure 3.3. 

Saunders exploits the form of these LU factors by maintaining a permutation 
of the columns and rows of U so that all the spikes appear on the right-hand side 
without changing their relative ordering, yielding: 



Figure 3.2: Spikes within a Bump 
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When the p th column of An is replaced by a k , the method proceeds as follows: 

1. Delete the p ,h column of U , u r , and move all the following columns one 
position to the left. 

2. Place o = L~ l a* in the right-hand column of U. 
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Figure 3.3: LU Decomposition of A n 

3. Move the p ,h row of {./, u p , to the bottom of {/, and move all rows in between 
up one position. Note that this preserves the triangular form of all but the 
last column of U, and further note that row u p has nonzeros only in columns 
corresponding to R. 

4. Using Gaussian elimination with row interchanges, transform u p to a sin- 
gleton column, thereby restoring the upper triangular form of U. Thus, 
U = E t . . . E\V , where U is the final updated form of U and E r . . . . ,E\ are the 
elementary transformation matrices corresponding to the Gaussian elimination 
steps (see, e.g.. Murtagh [1981]). 

5. Apply these same Gaussian elimination steps to L -1 , forming 

ir' = E t ...E x L~'. 



Typically, L~ l is held in product form and thus these transformations are 
simply added to the current list of transformation vectors (“eta” vectors). 

This method has the virtue that new nonzeros can be created only in the 
submatrix F of U, and thus R may be carried in peripheral storage. The numerical 
properties are reported to be quite good. 

The computational burden of continuously maintaining An in block triangular 
form is excessive, and the usual approach is to refresh the structure as part of the 
basis re-inversion routine. Between the re- inversions, the structure is left untended. 
An effective heuristic due to Hellerman and Rarick [1972] is commonly used for this 
purpose. 

A pivot in tableau coordinate ((n),(j)) results in a row exchange in An. 
Writing (3.6) as 

A u U~ l = L . (3.8) 

When row a T replaces the q th row of A u to form A n , we have 



A\iU~ i = 



~T _.j 

L \ 

X 



where row 0 — a T U 1 replaces the q th row of L. The desired structure of L may be 
restored by methods symmetric to those developed for the column exchange case, 
again forming : 



An = LU . 

A pivot in tableau coordinate ((H), (jj)) causes the dimension of A n to in- 
crease by one through the addition of a row and a column. It is convenient to add 
the new row to the bottom of A u and the new column to its right. If A u is the 
new kernel, then 



A u — 



An a k 

&r &r y k 



where a r ,a h and a r< k are original matrix coefficients. 

The desired updated decomposition is of the form: 



(3.9) 
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which requires that 

Lu k = a k 
l r U = a r , 



and 



&T,k — T 



k,k • 



Setting Uk,k = 1 • we have: 



lr,T = flr.it - Htl k 
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The final pivot coordinate, ((i), (j)) results in the dimension of A u decreasing 
by one. If An is the current transformation kernel and An results after the pivot, 
without loss of generality we have: 
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Q>r y k &r 




h,k 


0 




U r ,r 


u r 




1 

1 




l k 


Lx 




0 


U x . 



( 3 . 10 ) 



where (a r jt, a r ) is the leaving row, (a ri *, a k ) T the leaving column and a r> j t the 
pivot element. Using an analysis similar to Rosen [1960], we note that if the square 
matrix C is nonsingular and is partitioned as: 



C i C 2 
C3 C4 > 

where C\ is square and nonsingular, C 4 is square and C 0 = C 4 — CoCf 1 C 7 is non- 
singular, then 




C~ l = 



Cf 1 + C^C 2 Co X C 2 C ; 1 -C^C 7 Co l 
-Co'C’aCT 1 Cq 1 



Applying this result to (3.10), we discover that 



An o (u r ,fc) q t — L\U\ , 



or 



An — L\U\ + — J a k a T . (3.11) 

If the term o k a r should be the zero matrix, our new decomposition is at 

hand. L nfortunately. a k a r need not be zero, but we may guarantee it to 

be so by performing a preliminary column exchange update to A u ■ We replace 
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column ^ a r> k,d k ) with column (1,0)^ and compute the updated factors exactly as 
we did in the column exchange case considered earlier. This results in a modified 
transformation kernel A\\ with factors 



1 a T 




h,k 0 




U T ,r Ur 


i 

[r € 

o 




i 




0 Ui 



Then the second term of (3.11) is 




and thus the updated LU factors are at hand: 



An — L\tj\ . 

We now* see the advantages and disadvantages of the Graves mutual primal- 
dual method. It has the advantage that, although the transformation kernel can 
be as large as m, the number of structural constraints in the problem instance, its 
dimension is actually equal to the current number of binding structural constraints 
(or equivalently, to the number of potentially nonzero primal variables). At our 
initial basic solution, the kernel dimension is zero. If the maximum kernel dimen- 
sion during the course of the solution trajectory is much smaller than m, we enjoy 
significant computational advantages in storage, update and computation. For a 
great many LP model instances, this is precisely the case. 

The dynamic nature of the transformation kernel clearly complicates update 
procedures. While the usual primal Simplex method requires only the column ex- 
change update. Graves’ method requires four update cases. 
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In our implementation, we will seek methods for handling the transformation 
kernel that enjoy the advantages while mitigating the disadvantages. 

It is this general strategy for row and column generation we develop and en- 
hance in our specializations. By identifying a special structure in the problem data, 
we will introduce additional linear dependencies in the rows and columns of the 
tableau which will further simplify their generation and will also further reduce the 
dimension of the transformation kernel. 
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IV. FACTORIZATION 



A. INTRODUCTION 

Each of the algorithms we present in this research can be viewed as a special- 
ization of a general approach to large-scale linear programming developed by Graves 
and McBride [1976], which they call the “factorization approach”. It is based on 
distinguishing special rows and columns in a way that allows large parts of the ba- 
sic tableau to be represented implicitly, to be generated easily from the remaining 
explicit parts only when actually required by the algorithm. 

Although sometimes used interchangeably in the literature, we recognize a 
distinction between partitioning methods and factorization methods. In the former, 
a formal distinction is made between substructures which appear in the model in- 
stance, usually constraints or variables. An approach for solving the problem is then 
developed which exploits the substructures. For example, Dantzig- Wolfe decompo- 
sition is such an approach which partitions constraints, while in Benders decompo- 
sition the variables are partitioned. Such an approach may be applied statically, in 
which case the desired partitions are identified at the start of the solution process 
and remain fixed throughout, or it may be dynamic, in which case the partitioning 
may be adjusted as the solution proceeds. 

In contrast, we consider factorization methods to be those in which a formal 
distinction is made between substructures in the LP basis (and thus in the basic 
tableau). The algorithm is then specialized to exploit this special substructure. 
Thus, we view the GUB algorithm of Dantzig and Van Slyke [1967] and the mul- 
ticommodity network flow algorithm of Hartman and Lasdon [1970] as examples 
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of factorization algorithms. Factorization methods may also be static or dynamic, 
depending on whether the dimension of the factored substructures in the LP basis 
may vary during the solution process. 

Our research is concerned with dynamic row factorization algorithms. We 
develop the general dynamic factorization approach here and discuss the important 
aspects of the algorithm. In subsequent chapters, we specialize this general approach 
to each of several special LP row structures. This general development closely 
parallels that of Graves and McBride [1976]. 

B. THE FACTORED TABLEAU 

The problem to be considered is: 

(FLP) min : wy 

s.t. : Ey < r } explicit constraints 

Ey < b } factored constraints 
— Iy <0 } nonnegativity constraints , 

where y is a n vector of decision variables, w is a n vector of cost coefficients, E 
is an m by n matrix of constraint coefficients for “explicit” constraints with right- 
hand side m vector r, F is a p by n matrix of constraint coefficients for “factored” 
constraints with right-hand side p vector 6, and —I is the negative of the n by n 
identity matrix. In this general development, we refer to the E-type constraints 
as “factored" only to distinguish them from the “explicit” E-type constraints, and 
assume nothing about their structure. Not Until our specializations in later chapters 
will we impose special structure on E, and the structures w r e will consider permit 
the representation of the F - type constraints without the inversion of a matrix. We 
will see that this approach is centered around handling the part of the basis cor- 
responding to the E-type constraints explicitly while factoring the portion of the 



62 



basis corresponding to the .F-type constraints. The notation is chosen to suggest 
this idea. 

As we saw in Chapter 2, in the mutual primal-dual method the primal al- 
gorithm and the dual algorithm share the same tableau, and thus the tableau for 
(FLP) may be derived from either perspective. We present only the derivation from 
the primal viewpoint here. 

Recall that a basis for the primal algorithm consists of n linearly independent 
rows from the constraint coefficient matrix when it is assumed to include both 
structural (explicit and factored) and nonnegativity constraints. Assume that the 
current row basis consists of / rows from E, k rows from F and (n — k — l) rows 
from —I. Using the notation of Chapter 2 temporarily, we have: 



B = 



l+k^ 


n—l—k 






^ \ 




An 


A 12 


} / + k 


0 


-1 ) 


} n — l — 



where [An An] includes all basic structural rows, from both E and F. 

We will ultimately be interested in isolating the effect of each type of struc- 
tural constraint algebraically in the factored tableau, and thus we require greater 
resolution in our factored basis. Introducing obvious notation, we have: 



l 






[An A12] — 



En E\2 
Fn F 12 



n—l—k 

' \ 
Ei 3 

F 13 / 



}/ 

}k , 



where 



An] - 



En E\2 

Eh E\2 
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From the development in Chapter 2, we recognize An as the explicit kernel, 
and thus A n exists. It then follows that it is always possible to identify among the 
columns of [Fn Fi 2 ] a nonsingular submatrix Fu of dimension k, since otherwise 
the rank of [F n F 12 ] is at most (k — 1) and thus the rank of B is at most (n — 1). 
We will later see that one of the important implementation challenges is the task of 
efficiently managing the structure and nonsingularity of Fu- 
The full factored row basis is then: 



n-l-k 



B = 



Introducing the notation: 



U) 


£n 


£ 12 


£ 13 


}l 


in) 


£ n 


£12 


£13 


}k 


( jjj ) 


0 


0 


-I ) 


} n — l — 



( 4 . 1 ) 



■£'12 — -£11 Fjj 1 F 12 

£13 — e u f u 'f 13 

its inverse is: 




B- 1 



-Fu'/vAj (I + FjFuAjEnWn 1 F {? (£ 13 - F p A^ A n ) 

A l - l — FijFjj A l ^ A 12 

0 0 -I 



Grouping the coefficients of the nonbasic constraints and applying the same 
column ordering yields: 
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i — l — k 



D = 



(0 


-I 


0 


0 


}k 


(«) 


F n 


F 22 


F 23 


}p-k 


(Hi) 


0 


-/ 


0 


)l 


(iv) 


( P21 


D 22 


F23 ) 


} m — l 



We will explain the ordering of the rows of D shortly. 

The principal part of the factored tableau is DP, where P 
conjugate row basis. Introducing the additional notation: 



= -B 



-l 



F22 


A 


P 2 2 


- P 2 lP 11 1 Pl 2 


P23 


A 


P23 


-P 2 lP 1 l 1 Pl3 


^21 


A 


P22 


— P 2 iP n 1 Pi 2 


y4 2 2 


A 


E 23 


- P^Pr^Pia 



the principal part of the factored tableau is: 

(i) Uj) 



(0 

DP = <"> 
(Hi) 

(iv) 



— P 2 2 ^ 11 1 

i-i 
^11 

—AzxA^ 



{FnAh'En ~ F2i)Fn' 

—Au EuFu 1 
(Ai\A^i E\\ — E 2 \)F [ j 



(in) 



-Fu'FnA-J (I + FJ F n A-J E U )FJ 



Fn\F l2 -F u Ai{A 

P?3 — P22^4]l 1 ^12 



atU 



11 ^12 



^22 — A 2 \A^ A \2 



Partitioning w = (u'i . w 2 . w 3 , w 4 ), r T = (r 1 ,r 2 ) T and b T = (6j, 6 2 ) 
introducing the notation: 



U'2 



= u> 2 - w i F n 1 F n 



( 4 . 2 ) 



is the 



\ 

12) 



/ 

( 4 . 3 ) 

T and 
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U’3 = U) 3 -W 2 F u 1 Fi3 

b 2 = bt-FnFjbx 
r, = r.-EuF-H, 
f 2 = r 2 , 



the complete factored tableau is: 



-Fi/Fn/t,, 1 

— Fj2y4ii’ 

i-i 

. n . 

— /An/if, 1 
-wj/in 1 



(/ + Fn'FnAXEn)!# 

j 1 Fn — FjjJF,, 1 
— /in'FnF,! 1 
(/Ci/tjYFn — Eu)F n ' 

(tl'j/t ji'Fji — W\)F id 



FfJ *_( Fa ~ F n A 7/ A l2 ) 

F 2 3 —-F 12 Aii A 12 

A22 — /4 2 iy4J' 1 1 i4i2 
tZ> 3 - wjA^Au 



F n ](bi - FuAi'ri) 
62 — F22>4f, 1 f l 

ir.'r, 

fj - 4,1 ir. 1 /. 

U’lFi’&i + 



(4-4) 



We see from (4.4) that with knowledge of the current factorization, we can 
construct the entire tableau from F,^ 1 , A f/ and the original problem data. The 
dimension of Ffj 1 is equal to the number of F-type constraints that are currently 
basic, and thus can be at most p. The dimension of /If/ is equal to the number of 
F-type constraints that are currently basic. Hence, its dimension cannot exceed m. 
We call Af/ the explicit transformation kernel and F^ 1 the factored transformation 
kernel. 

We have seen in Chapter 3 that the origin is always a convenient initial basic 
solution for the mutual priinal-dual method, and the same is true for the factoriza- 
tion approach. An initial basic solution is always: 
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B° = [-/] 




Starting from this solution, we may view the algorithm as progressing by exchanging 
rows of F and E from D with nonnegativity constraints from B. We will find it useful 
to associate with each F-type constraint in B a unique nonnegativity constraint in 
D. Borrowing the terminology of Dantzig and Van Slyke [1967], we will call each such 
unique nonnegativity constraint (and its corresponding variable) the “key” variable 
for the associated F-type constraint. The algebraic structure arising from this basic 
F-type constraint /nonbasic “key” variable association allows us to represent parts 
of the tableau implicitly rather that explicitly. To distinguish these “key” variables 
in D from others, we place them in region (?) of D. This is the reason for the row 
ordering mentioned earlier. 

C. BENEFITS OF FACTORIZATION 

Graves and McBride [1976] report three principal benefits of the factorization 
approach: 

1. Good starting bases; 

2. Reduction in memory requirements; 

3. Reduction in work per pivot. 

Although the origin is always available as an initial basic solution, Graves and 
McBride [1976] suggest that a better one can frequently be found with only small 
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computational effort using this approach. They form the Lagrangian Dual of (FLP) 
with respect to the E - type constraints, we obtain: 

(LFLP) min : wy + A(Ft/ — r) 

s.t. : Fy < b 

-Iy< 0 

where A > 0 is a m vector of dual variables. If A equals A*, the optimal dual variables 
for (FLP), then an optimal solution of (FLP) is among the optimal solutions of 
(LFLP). They speculate that if A is a “good” estimate of A*, then using A in (LFLP) 
and solving for y should yield a “good” starting point y for (FLP). If no such estimate 
A is available, then A = 0 may be used. Depending upon the structure of the F-type 
constraints, the resulting problem may either be much easier to solve than (FLP) and 
thus may be solved optimally, or heuristics may be used to find a computationally 
inexpensive “good” solution. All factorization algorithms can be started by solving 
(LFLP) first. 

Standard simplex techniques require a basis of dimension (m + p ) to solve 
(FLP). The factorization approach requires the explicit transformation kernel Ay j 1 , 
whose dimension is at most m, and the factored kernel Fu , whose dimension is 
at most p. Thus, it is clear that we expect a considerable reduction in memory 
requirements using the factorization approach. When we specialize the approach 
for models in which the F-type constraints have special structure (e.g., GUB, pure 
network or generalized network), we will see that memory requirements can be 
further reduced. 

While it is difficult to measure the computational effort per pivot when product 
form inverse or basis decomposition techniques are used, there are indications that 
the factorization approach can yield substantial benefits. Both analytic (McBride 
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[1972]) and empirical (Tomlin [1972]) studies indicate that when the F- type con- 
straints are GUB constraints (we will cover this in detail in Chapter 6), the factoriza- 
tion approach results in per-pivot computational benefits when the GUB constraints 
comprise at least 30% of the total number of structural constraints. With some sim- 
plifying assumptions, a similar analysis can be developed for the case where the 
F-type constraints are pure network rows. 

We will consider each of three different factored row structures in this work. 
The complexity of the factored row structure determines the computational difficulty 
of representing and updating the factored kernel F u . The simplest row structure 
we consider, in which the factored rows are generalized upper bounds, allows F n to 
be interpreted as a simple permutation matrix. The second factorization, in which 
the F-type constraints are pure network rows, is a more complex structure which 
subsumes generalized upper bounds as a special case. In this case, Fn may be 
interpreted as representing a partial ordering defined over a subset of the rows of F. 
The final factorization, in which the rows of F are generalized network rows, is the 
most general we consider in this w'ork and subsumes both the other factorizations. 
Here, Fn may be interpreted as a linear transformation in which a near partial 
ordering exists among a subset of the rows of F. 
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V. GENERAL IMPLEMENTATION TOOLS 



A. INTRODUCTION 

We now provide an overview of our implementation of the factorization method. 
We describe the representation and update of the three important algorithm compo- 
nents: the basic tableau, the factored kernel and the explicit transformation kernel. 
The fundamental update steps are described in terms of functions which operate 
on each of the algorithm components. Where the details of the update actions are 
common to all three factorizations treated in this research, the details are presented 
in this chapter; otherwise, they appear in subsequent chapters. 

B. THE FACTORED TABLEAU 

As before, we are interested in the problem: 



min : wy 

y 

s.t. : Ey < r } explicit constraints 

Fy < b] factored constraints 
—Iy < 0 } nonnegativity constraints 



Recall the algebraic representation of an arbitrary primal row basis: 
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E n 


E 12 


E 13 
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-I) 
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with the corresponding nonbasic rows: 



\ 



(0 


-1 


0 


0 


(») 


Fn 


£ 22 


£23 


(*«) 


0 


-/ 


0 


(iv) 


v £21 


£ 22 


£23 / 



The conjugate row basis inverse is: 



(5.2) 



(5.3) 



P = -B _1 = 



l FnFvA’J 
_ 2- 1 

0 



-Fj -F^FvAu'EnFj 
A\i FiiFn 1 
0 



-Fii'Fn + FZ'FtiArflEn- EnFtfFn) 

— /4 1 1' ( £^i 3 — E\\F n l F\z) 

1 



and the principal part of the factored tableau is: 
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-Fn'FaA^' 


FiV + F7i l FnA X i EuF u 
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— F?iA xx + 


FnAii E U F 'ii* - F- 2 \F u l 






-FuFTSFuAu'EnF^ 
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(5.4) 



Ff.'F.j - Fu'F l} A^'(E l3 - EnFtfFn) 

F 73 - FjiF, 1 'F ,3 + Fj.Fri'Fn^n'fFn - ^..Ff.'F.j) 
-F n A~'(E l3 - Fuff/ F n) 
ir.'(Fi 3 -F..F,': , F 1 3) 

F 7 j F E 7 i F’ i| 1 / r 'j2^|| (£*13 FiiFu F13) 

-FjiFf.'F.j - £»*»(£» - FnFf.'Fn) ; 
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The entire principal part of the factored tableau can be constructed with 
knowledge of the row and column ordering, Fn (or F n ) and . Our discussion 
of the implementation centers around the representation and update of these three 
components. 

We restrict our attention to the principal part of the tableau because our 
implementation maintains a current representation of the right-hand side column 
and bottom (cost) row (collectively referred to as the “problem rim”) at all times. 
These quantities are always immediately available and need not be computed upon 
demand. This is in contrast to (5.4), where only and F\\ are always available 
and any other quantity of interest (for example, a complete row or column) must 
be computed. 

C. ALGORITHM OVERVIEW 

In broad terms, the algorithm proceeds as follows: 

1. Establish an initial basic solution. 

2. Stop if the current solution is optimal. Optimality exists when the right-hand 
side column is nonnegative and the bottom (cost) row is nonpositive. 

3. Based on a pricing scheme and a ratio test, select and generate a row and 
column of (5.4). 

4. Stop if this row and column reveal infeasibility or unboundedness. 

5. Transform the problem rim. the tableau representation, F u and A^ 1 , and 
return to step 2. 

After developing the data structures and necessary update operations, we will 
expand upon this overview. 
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D. IMPLEMENTATION CONVENTIONS 



We will first discuss the implementation of each of the three algorithm com- 
ponents individually, and then provide a detailed presentation of the complete al- 
gorithm. Our implementation is coded in the FORTRAN language, and thus the 
following discussion is colored by the character of that language. 

1. The Tableau 

The mutual primal-dual method fosters a symmetric view of the linear 
programming problem and the solution process. This symmetry is reflected in the 
indexing conventions required in our implementation. Assuming now that (5.1) 
consists of m structural constraints (both factored and explicit) and n structural 
variables, we require the structural constraints to be indexed from 1 to m, and the 
variables from (m + 1) to (m -f n). While our statement of the problem specifies 
the nonnegativity constraints explicitly, we will of course deal with them implicitly. 
Note, however, that we have two equivalent interpretations of the indices (m -f 
1) through (m -f n). We can view them as the indices of the structural problem 
variables or as the indices of the nonnegativity constraints corresponding to each 
of the structural variables. We will use both interpretations in our development. 
Observe also that we require no explicit representation of logical (slack, surplus or 
artificial) variables. 

The indexing of the original problem data exists independently of the 
solution process and is thus referred to as the “extrinsic” indexing. The factoring of 
the constraints to form B and D represents a second indexing of the problem which 
defines the current solution, independent of the extrinsic indexing. We refer to this 
indexing system as “intrinsic”, and it partially defines the current representation of 
B and D. Our convention for intrinsic indexing has the rows of D labeled from 1 
to m and the rows of B labeled from (m + 1) to (m + n). 
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The original problem data is stored in super-sparse form, (in super-sparse 
form, each real value is stored once, and each nonzero coefficient is represented by a 
row index, a column index, and a pointer to a real value) with unique nonzero real 
values stored once in an array of type DOUBLE PRECISION, and referenced by the 
FORTRAN equivalent of a pointer from each constraint matrix coefficient. These 
pointers to the real values are singly-linked by both row (with associated column 
indices) and by column (with associated row indices). This allows symmetric row- 
wise and column-wise access to the data to the performed efficiently, which is an 
important virtue in a primal-dual design. 

To specify the structure of (5.4) we require a representation of the cur- 
rent mapping between the extrinsic indexing of the problem data and the intrinsic 
indexing of the current tableau. We represent this mapping by two arrays of type 
INTEGER, each dimensioned from 1 to (m + n), which are used as pointers. We 
denote these arrays MINT() and MEXT(). MINT() represents the mapping from 
intrinsic to extrinsic indices, and MEXTQ represents the inverse mapping. Thus, 
if IINT is an intrinsic index of a row or column of (5.4) and IEXT is the extrinsic 
index of the constraint or variable currently mapped to position IINT of the tableau, 
then: 



MINT(IINT) 


= IEXT 


MEXT(IEXT) 


= IINT 


MINT(MEXT(IEXT)) 


= IEXT 


MEXT(MINT(IINT)) 


= IINT. 



To complete the representation of (5.4) we require knowledge of the fac- 
torization which is represented algebraically by regions (i). (n). (Hi), etc. We 
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handle this with the use of several variables of type INTEGER which are used as 
pointers to record the ending intrinsic index of the various tableau regions. Table 
(5.1) lists the regions and the associated pointer names. 

TABLE 5.1: Indices of Tableau Regions 
REGION BEGINNING INDEX ENDING INDEX 



0) 


1 


MKC 


(ii) 


MKC 4- 1 


MFR 


(iii) 


MFR + 1 


MXR 


(iv) 


MXR + 1 


M 


(j) 


M + 1 


NXR 


(jj) 


NXR + 1 


NFR 


(jjj) 


NFR + 1 


M + N 



The arrays MINTQ and MEXT() and the pointers MKC, MFR, MXR, 
NXR and NFR comprise the data structures necessary to represent the principal 
part of the tableau. To complete the representation of the tableau we require a rep- 
resentation of the right-hand side column, the bottom row and the current value of 
the extremal function. The right-hand side column and bottom row are represented 
by an array of type DOUBLE PRECISION, dimensioned from 1 to (m-|-n), referred 
to as RIM(). Positions 1 through m correspond to the right-hand side column and 
positions (m + 1) through (m . + m) correspond to the bottom row. Note that these 
are intrinsic indices. The current value of the extremal function is held in a scalar 
variable of type DOUBLE PRECISION called OPT. 

We now discuss the operations necessary to maintain and update the 
tableau representation. To recognize an optimal solution we simply scan the RIM() 
array. When the values in positions 1 through m are nonnegative and those in 
positions (m + 1) through (m + n ) are nonpositive, the current solution is optimal. 



10 



Violations of either sign discipline may serve as potential primal algorithm target 
rows or dual algorithm target columns, respectively. 

Having selected a target row or column, we must generate the correspond- 
ing row or column of the tableau, compute ratios with respect to the right-hand side 
column or bottom row, and thus identify the index of a corresponding column or 
row. We then generate that column or row of the tableau. Because the details of row 
and column generation differ according to the factorization, we defer their discussion 
until the appropriate chapter. To represent the update operations, however, we de- 
fine the abstract functions Generate_Row(IPRI) and Generate_Column(IPCI). 
These functions accept as arguments the intrinsic row (IPRI) or column (IPCI) 
index of the current tableau and return the current representation of that row or 
column. We require a working array of type DOUBLE PRECISION to hold these 
current representations. We define ROWCOLQ to be such an array, dimensioned 
from 1 to (m + n), where positions 1 through m correspond to column IPCI of (5.4) 
and positions (m + 1) through (m + n) correspond to row IPRI. ROWCOL0 is 
indexed intrinsically in conformance with our convention for (5.4). 

Recall that our interpretation of a pivot from the primal algorithm view- 
point is as an exchange of rows between B and D. The row indices of D correspond 
to the row indices of (5.4), while the row indices of B correspond to column indices 
of (5.4). Thus, a pivot requiring the exchange of a row of B and a row of D requires 
the exchange of a row index of (5.4) with a column index of (5.4). Using MINT() 
and MEXT(), such an exchange requires just a few assignment statements. For 
example, if the pivot row IPRI corresponds to extrinsic index EPRI and the pivot 
column IPCI corresponds to extrinsic index EPCI, then prior to a pivot we have: 
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MINT(IPRI) = EPRI 
MEXT(EPRI) = IPRI 
MINT(IPCI) = EPCI 

MEXT(EPCI) =IPCI 
and after the pivotal exchange we have: 

MINT(IPRI) = EPCI 
MEXT(EPCI) = IPRI 
MINT(IPCI) = EPRI 
MEXT(EPRI) =IPCI. 

We call this type of index exchange a “primary exchange”. For abstrac- 
tion purposes we define a function Primary _Index_Exchange (IPRI,IPCI) which, 
given the current tableau representation and the intrinsic pivot row index IPRI and 
the intrinsic pivot column index IPCI, performs the appropriate update of MINTQ 
and MEXT(). Every pivot requires a primary exchange. 

In addition to the partitioning of constraints into B and D , we are re- 
quired to maintain the factorizations ((i).(ii), etc.) within each. Having performed 
the primary exchange, it is possible that additional exchanges within B and/or D 
are required to restore the proper factorization. For example, if the primary ex- 
change removes an F-type constraint from D (region («)) and places it in region ( j ) 
of B , an additional exchange is necessary to move the row from region (j) to region 
(jj ) . We refer to such actions as “secondary exchanges”. Secondary exchanges in D 
correspond to row exchanges in (5.4) and secondary exchanges in B to column ex- 
changes in (5.4). Thus, we define the functions Row’ Jndex _Exchange(IRIl .IRI2) 



and Column _Index_Exchange(ICIl,ICI2) which exchange intrinsic rows IRIl and 
IRI2 or intrinsic columns ICI1 and ICI2, respectively. Since the indexing of the 
RIM() array corresponds to the tableau indexing, we assume that any secondary 
exchanges performed on the tableau are also performed on RIM(). 

Finally, our factorization requires that the factored kernel, Fn, remain 
nonsingular. It is possible that the primary and secondary exchanges will destroy 
the nonsingularity of Fn. When this occurs, additional row exchanges between 
regions (i) and (m) of D will be necessary to restore the required nonsingularity of 
F\\. We call such exchanges “tertiary exchanges”. Again, we assume that tertiary 
exchanges performed on the tableau are also performed on RIM(). 

The factorizations of B and D are dynamic in nature, meaning that 
a particular pivot may cause the dimension of one or more regions of B and/or 
D to increase or decrease by one row. Such dimension changes are handled by 
simply adjusting the values of the factorization pointers. We define the functions 
Increment (XXX) and Decrement(XXX), where “XXX” is MKC, MFR, MXR, 
NXR or NFR, to increment or decrement, respectively, the appropriate pointer 
value. 

In summary, the tableau data structures are: 

MINT() 

MEXT() 

RIM() 

ROWCOL() 

MKC. MFR, MXR, NFR, NXR, 
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and the necessary update operations are: 



Generate JRov/ (I PRI) 

Generate_Column(/PC7) 

Primary Jndex_Exchange(//\ft/, I PCI) 
Row_Index_Exchange(//?/l, IRI2) 

Column Jndex_Exchange(/C71,/C72) 

Increment(A r A r A r ), and 
Decrement(JVAW). 

2. The Explicit Transformation Kernel 

As can be seen from (5.4), the rows of A f, 1 correspond to region (Hi) of 
the tableau (nonbasic nonkey variables) and the columns to region (j) (basic explicit 
constraints). The dimension of AJ"/ may vary as the algorithm progresses, and the 
mechanism for these changes is the addition and deletion of rows and columns. 
We thus require data structures that permit the efficient insertion, deletion and 
exchange of rows and columns of A^ 1 . Further, the number of nonzero elements 
in A 7 / varies from pivot to pivot, independently of dimension changes. We store 
these nonzero elements in a stack, referring to them via pointers which are stored 
as an orthogonally linked list, doubly linked by row and doubly linked by column. 
This allows the efficient update of Af, 1 and accommodates the dynamics, at some 
expense in storage. Contrary to the tableau representation, we maintain the explicit 
transformation kernel representation indexed extrinsically rather than intrinsically. 
This allows permutations within region ( j ) and/or within region (iii) of the tableau 
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without requiring corresponding permutations in the columns and/or rows of the 
AJ representation. 

We again define a suite of functions which operate on our representation of 
A^ 1 and describe its update. Add_Row(IRE) and Add_Column(IJE) append ex- 
trinsic row IRE and extrinsic column IJE, respectively, to our representation of A^ 1 . 
We will ensure that the current representation of each is in work array ROWCOL() 
prior to executing these functions. Columns will always be appended to the right- 
hand side of Aii an d rows will be appended to the bottom. Delete_Row(IRE) and 
Delete_CoIumn(IJE) delete extrinsic row IRE and IJE, respectively, from the rep- 
resentation of Ail . Replace_Row(IREl,IRE2) causes the overlay of extrinsic row 
IRE1 by extrinsic row IRE2. The representation of IRE2 will have been previously 
placed in ROWCOLQ. Finally, Update_Explicit_TVansformation_Kernel per- 
forms a pivot update of the representation of Aj"/, using the current representation 
of the pivot row and pivot column in ROWCOL(). 

3. The Factored Kernel 

The data structures and update actions for the factored kernel vary ac- 
cording to the factorization and their discussion will be deferred. Here we define the 
necessary abstract update functions. 

Is_Factored_KerneLSingular tests the current representation of F u 
for singularity and returns the appropriate Boolean value. It will be used in certain 
pivot cases to determine which of several courses of action should be taken. 

All tertiary exchanges take the form of row exchanges between regions 
(?) and (???) of (5.4). Note that region (?) consists of the nonnegativity constraints 
of key variables. By our alternate interpretation of the nonnegativity constraint 
indices we may consider these to be key variable indices. Since the columns of F n 
are key variables, we interpret region (?) as containing the indices of the columns 
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of Fix* An exchange of rows between (t) and (in) of (5.4) is then interpreted as an 
exchange of columns of *h. 

Before such a column exchange for F n (and thus a row exchange in (5.4)) 
can be made, we must frequently identify the index of the column to be removed from 
Fn (which is an index in region ( i ) of (5.4)) and the index of the column which will 
replace it (which is an index in region (in) of (5.4)). We thus define Find_Column 
_to_Remove(IOUT), which selects from among the indices of region (i) of (5.4) 
the intrinsic index IOUT of the column to be removed from Fn. Similarly, Find 
_Column_to_Add(IIN) selects from among the indices of region (Hi) of (5.4) the 
intrinsic index IIN of the column to be be added to Fn- Finally, Update JFactored 
.Kernel updates the data structures used to represent Fn. 

E. COMPLETE ALGORITHM DESCRIPTION 

We are now in a position to fully describe our implementation of the 
algorithm. We do so by expanding the discussion of each step given previously in 
section C. 

1. Initialize the algorithm using the origin as the initial basic solution. Then 



B° = (-/) 




2. Stop if the current solution is optimal. The current solution is optimal if 
(RJM(IR) > 0 for 1 < IR < M and RIM(JC) < 0 for M+l < JC < M+K. 

3. (a) Select a violated primal or dual constraint as the target row or column. 

respectively. For purposes of exposition, assume the current solution is 
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primal feasible and we are executing the primal algorithm. Then the 
target row is the bottom row, and we select the intrinsic index IPCI of a 
column which will allow us to make a gain in the value of the extremal 
function (i.e., RIM (I PC I) > 0). 

(b) Generate_Column(IPCI). 

(c) Calculate ratios by computing RIM(IR)/ROWCOL(IR ) for those 1 < 
IR < M with ROWCOL(IR) > 0. UROWCOL(IR) < 0 for all IR, 1 < 

IR < M stop, since the problem is primal unbounded. Otherwise, select 

IPRI to be the intrinsic index yielding the smallest such ratio. Break 
ties in accordance with the following hierarchy: region (in) first, then 
region (n), then (?) and finally region (iv). Within a region, break ties 
by selecting IPRI to be the first such index encountered. 

(d) Generate_Row(IPRI). 

4. If, contrary to our assumption in 2. a, the current solution is not primal feasible, 
the target row would be some row IPRI of (5.4), 1 < IPRI < M , rather than 
the bottom row. We would proceed by next executing Generate_Row(IPRI). 
If ROWCOL(JC) > 0 for all M -f 1 < JC < M + N, we would conclude that 
the problem is primal infeasible. 

5. (a) Primary _Lnd ex .Exchange (IPRI, IPCI). 

(b) Pivot update RIM() and OPT using ROWCOL(). 

(c) Perform the necessary secondary and tertiary exchanges, as shown in Ta- 
ble 2. Note that some tertiary update actions are conditional, depending 
upon the singularity of F n . We use a notation similar to the FORTRAN 
BLOCK IF statement to indicate the conditional actions. 
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(d) Update_Factored_Kernel. 

(e) Update_Explicit_Transformation -Kernel. 

(f) Go to Step 2. 

We now give a more detailed explanation of the secondary and tertiary 
exchanges listed in Table 5.2, discussing each pivot coordinate individually. 

1. Pivotal coordinate ((i), ( j )). EPCI, the extrinsic index of an E-type constraint 
located in position IPCI of B, is exchanged with EPRI, the extrinsic index 
of a nonnegativity constraint located in position IPRI of D. The initial pivot 
exchange places EPCI in position IPRI of ( i ) and EPRI in position IPCI of 
(j), and thus secondary row and column exchanges are necessary to restore the 
structure of the tableau. Beginning with the column structure, we exchange 
column EPRI in position IPCI of region (j) with the extrinsic index located in 
position NXR of region ( jjj ), and then decrement the value of NXR. This has 
the effect of shifting the starting column index of region (jj) one position to 
the left. Now extrinsic index EPRI is in position (NXR-f-1) of the tableau and 
thus is still misplaced. We then exchange EPRI in position (NXR+1) with the 
extrinsic index located in position NFR. This places EPRI contiguously with 
region (jjj), so by decrementing NFR we restore the column factorization. 
The effect has been to decrease the dimension of (j) by one, increase the 
dimension of (jjj) by one and to shift the entire (jj) region one position 
to the left. In the row structure, EPCI in position IPRI of region (i) must 
eventually be moved to region (iv). However, prior to this pivot EPRI was a 
key variable, located in region (i). Its removal has disturbed the structure and 
thus the nonsingularity of F n . We must therefore identify a column currently 
in region (iii) which can be used to restore the nonsingularity of F n . We 
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thus invoke Find_Column_to_Add(IIN), which identifies the intrinsic index 
UN of such a column. Suppose the extrinsic index of the column in position 
IIN is KX. We exchange EPCI in position IPRI of region (i) with KX in 
position IIN of region (Hi). We now have extrinsic index EPCI in position 
IIN of (in') and we must move it to region ( iv ). We thus exchange EPCI 
in position IIN in region (in) with the extrinsic index in position NXR of 
(Hi), and then decrement NXR. This completes the restructuring of the row 
factorization with the effect of decreasing the dimension of region (Hi) while 
increasing that of region (it>). The effect of these exchanges has been to reduce 
the dimension of A fj 1 through the deletion of column EPCI and row KX. To 
maintain the proper structure of , we therefore Delete_Row(KX) and 
Delete_Column(EPCI). 

2. Pivotal coordinate ((z), (jj)). EPCI, the extrinsic index of a basic F-type con- 
straint located in position IPCI of B , is exchanged with EPRI, the extrinsic in- 
dex of a nonbasic nonnegativity constraint located in position IPRI of D. Since 
both EPCI and EPRI are misplaced after the primary exchange, secondary row 
and column exchanges are needed. In the columns, we exchange EPRI in po- 
sition IPCI of region (jj) with the extrinsic index located in position NFR of 
region (jj). Decrementing NFR then completes the column exchanges, leaving 
the dimension of region (jj) reduced by one and that of (jjj) increased by 
one. In the rows, we exchange EPCI in position IPRI of region (i) with the ex- 
trinsic index in position MI\C of region (z). Decrementing MKC then restores 
the row factored structure of the tableau. It remains to be seen, however, 
whether what remains of F n after the removal of row EPCI and column EPRI 
is still nonsingular. We test F n by invoking Is_Factored_Kernel_Singular. 
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If the response is FALSE, then no tertiary exchanges are required. If, how- 
ever, the response is TRUE, we must identify a column exchange that we will 
restore the nonsingularity of Fu. To identify the two columns which will be in- 
volved in the exchange, we invoke Find_Column_to_Remove(IOUT) which 
returns the intrinsic index IOUT (with associated extrinsic index NX, say) of 
a column to be removed from F u . We then call Find_Column_to_A.dd(IIN) 
which identifies the intrinsic index IIN (with associated extrinsic index KX, 
say) of a column in region (Hi) which can be used to replace IOUT. We now in- 
voke Row _IndexJExchange(IIN, IOUT), which completes the update of the 
row factorization. This row exchange results in the replacement of row KX of 
by row NX of the tableau. We currently have no representation of row KX of 
the tableau, so we compute it by invoking Generate_Row(IOUT). We may 
then restore the structure of Ajj 1 by invoking Replace_Row(KX,NX). 

3. Pivotal coordinate ((i), (jjj))- EPCI, the extrinsic index of a basic nonneg- 
ativity constraint located in position IPCI of B, is exchanged with extrinsic 
index EPRI located in position IPRI of D. The primary exchange properly 
places EPRI in region (jjj), and thus no secondary column exchanges are 
needed. Likewise, EPCI is properly placed in position IPRI of region (t), so 
no secondary row exchanges are required. Prior to the pivot, EPRI was desig- 
nated a key column and its removal disturbed the structure of Fji . Since EPCI 
is itself a column, it is possible that it may be a replacement for EPRI. To 
determine if this is the case, we invoke Is_Factored_Kernel_Singular. If the 
response is FALSE, the direct exchange of column EPCI for column EPRI in 
is justified and the tableau update is complete. If the response is TRUE, how- 
ever. a tertiary row exchange is required. EPCI, currently located in position 
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IPRI of region (z), will be the column we remove from F n , and its replacement 
is found by invoking Find_Column_to_Add(IIN) which returns the intrinsic 
index IIN (associated with extrinsic index KX, say) located in region (iii). We 
perform a tertiary exchange by invoking Row _IndexJExchange(IIN, IPRI). 
This completes the row factorization update. This tertiary exchange results 
in the replacement of row KX of A ^ by row EPCI of the tableau. We must 
therefore perform the corresponding row exchange in A^ also. To perform 
such an exchange, we require the current representation of row IPRI of the 
tableau. This is precisely the pivot row and is already available in the work 
array ROWCOLQ. Thus, we may invoke ReplaceJRow(KX,EPCI) directly. 

4. Pivotal coordinate ((H), (j))- EPCI, the extrinsic index of a basic E-type 
constraint located in position IPCI of B , is exchanged with EPRI, the ex- 
trinsic index of a nonbasic F-type constraint located in position IPRI of D. 
After the primary exchange, both EPCI (in (z'z)) and EPRI (in (j)) are mis- 
placed, and thus both secondary row and column exchanges are necessary. 
We exchange EPRI in position IPCI of region (j) with the extrinsic index 
in position NXR of region (j) and then decrement NXR. This has the ef- 
fect of reducing the dimension of region (j) while increasing that of region 
(jj). Notice that we have added row EPRI to the structure of E ll5 and 
thus an additional column must at some point be identified. Several row 
exchanges are needed. First, we exchange row EPCI in position IPRI of 
region (ii) with the extrinsic index located in position (MKC+1) of region 
(z’z) and then decrement MKC. This restores the structure of region (ii) by 
reducing its dimension by one row. We now must execute a series of ex- 
changes that ultimately places EPCI in region (iv) and also adds a column 
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to region ( i ) to restore the nonsingularity of F n . To do this, we first iden- 
tify the column to be added to F n by invoking Find_Column_to_Add(IIN), 
which identifies the intrinsic index IIN in region (n’t) (associated with an 
extrinsic index KX, say) which may be added to Fn. Notice that the dimen- 
sion of Fn is increased by one through the addition of row EPRI and col- 
umn KX. We then perform Row Jndex_Exchange(IIN,MKC), which places 
KX in position MKC of region (t) (correctly) and places EPCI in position 
IIN of region (Hi) (incorrectly). To complete the tertiary row exchanges, we 
invoke Row _Index_Exchange(IIN,MXR), followed by Decrement(MXR). 
This properly places EPCI in region ( iv ) by increasing the dimension of re- 
gion (iv) by one row while decreasing that of region (in). Finally, notice that 
the dimension of An has been reduced by one through the removal of column 
EPCI and row KX. To update the representation of An accordingly, we invoke 
Delete_Row(KX) followed by Delete_Column(EPCI). 

5. Pivotal coordinate ((ii), (jj))- EPCI, the extrinsic index of a basic F-type 
constraint located in position IPCI of B , is exchanged with EPRI, the ex- 
trinsic index of a nonbasic F-type constraint located in position IPRI of D. 
The primary exchange places both EPRI and EPCI in their proper positions, 
so no secondary row or column exchanges are needed. However, the struc- 
ture of has been altered through the addition of row EPRI and the deletion 
of row EPCI. To determine if the structure and nonsingularity of F n has re- 
mained intact, we invoke Is_Factored_Kernel_Singular. If the response is 
FALSE, the tableau update is complete. If the response is TRUE, we must 
identify a column exchange that will restore nonsingularity. Thus, we invoke 
Find_Column_to_Remove(IOUT) which returns the intrinsic index IOUT 
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(associated with the extrinsic index NX, say) of a column located in region 
(i) which may be removed from Fn- We invoke Find_Column_to_A.dd(IIN) 
which returns the intrinsic index (associated with an extrinsic index KX, say) 
of the column IIN located in region (in) which may be added to Fn. We then 
perform a tertiary exchange by invoking Row_Index_Exchange(IIN,IOUT). 
This last action results in the replacement of row IIN of A n* by row IOUT of 
the tableau. We have no current representation of row IOUT, and to compute 
one we invoke Generate_Row(IOUT). We then update the representation of 
A n by invoking Replace_Row(KX,NX). 

6. Pivotal coordinate ((H), (jjj))- EPCI, the extrinsic index of a basic nonnega- 
tivity constraint located in position IPCI of B , is exchanged with EPRI, the 
extrinsic index of a nonbasic F-type constraint located in position IPRI of D. 
The primary exchange misplaces both EPRI and EPCI and thus both sec- 
ondary row and column exchanges are necessary. We exchange EPRI, located 
in position IPCI of region (jjj) with the extrinsic index located in position 
NFR. By then incrementing NFR, we restore the column factorization struc- 
ture by increasing the dimension of region (jj) while reducing that of region 
(jjj). Notice that the addition of row EPRI to the structure of Fn indicates 
that an a corresponding column must also be located. To restore the structure 
of region (ii), we exchange EPCI in position IPRI of region ( ii ) with the extrin- 
sic index in position MXC of region (i). Incrementing MKC restores the struc- 
ture of region (ii) by reducing its dimension by one row. Since EPCI is a col- 
umn which we have now placed in position MKC, it is possible that Fn is now 
nonsingular. To determine this, we invoke Is_Factored_Kernel_Singular. 
If the response is FALSE, the tableau is complete. If, on the other hand, 
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the response is TRUE, we must perform a tertiary column exchange to re- 
store the nonsingularity of F u . We will be removing EPCI, currently located 
in position MKC, from Fij and replacing it with the column identified by 
Find_Column_to_Add(IIN), which is extrinsic index KX located in posi- 
tion IIN of region ( iii ). We then perform Row _Index_Exchange(IIN,MKC), 
which completes the tertiary exchange. This exchange implies a row exchange 
in Aii also. We are required to replace row KX of Aii with the current 
representation of row EPCI of the tableau. Since IPRI is the pivot row, the 
current representation of EPCI is readily available in ROWCOL0. We may 
then complete the row exchange by invoking Replace_Row(KX,EPCl). 

7. Pivotal coordinate ((m),(j)). EPCI, the extrinsic index of a basic E-type 
constraint located in position IPCI of B, is exchanged with EPRI, the extrin- 
sic index of a nonbasic nonnegativity constraint located in position IPRI of 
D. The primary exchange misplaces both EPRI and EPCI, and thus both 
secondary row and column exchanges are necessary. We exchange EPRI, in 
position IPCI of region (j), with the extrinsic index in position NXR of re- 
gion (j). We then exchange EPRI, now in position NXR of region ( j ), with 
the extrinsic index located in position NFR of region (jj). We complete the 
column update by decrementing both NXR and NFR, which has the effect of 
reducing the dimension of region (j) by one column, shifting region (jj) one 
column to the left in the tableau and increasing the dimension of region (jjj) 
by one column. For the row update, we first exchange EPCI, in position IPRI 
of region (iii), with the extrinsic index located in position MXR of region 
(iii). By then decrementing MXR we restore the row factorization by increas- 
ing the dimension of region (iv) by one row while decreasing that of region 
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(tit). The effect on A^l has been to reduce its dimension by one through the 
deletion of column EPCI and row EPRI. The corresponding update to the 
representation of A ^ requires that we invoke Delete_Row(EPRI) followed 
by Delete_Column(EPCI). 

8. Pivotal coordinate {(in), (jj))- EPCI, the extrinsic index of a basic F-type 
constraint located in position IPCI of B , is exchanged with EPRI, the extrin- 
sic index of a nonbasic nonnegativity constraint located in position IPRI of 
D. The primary exchange misplaces both EPRI and EPCI, and thus both 
secondary row and column exchanges are necessary. We exchange EPRI, cur- 
rently located in position IPCI of region (jj), with the extrinsic index currently 
located in position NFR of region (jj). By then decrementing NFR, we re- 
store the column factorization by decreasing the dimension of region (jj) while 
increasing that of region (jjj). Note that the removal of EPCI from the struc- 
ture of F n implies that the dimension of F n decreases by one, and thus a 
corresponding column must be removed from F u ■ To locate this column, we 
invoke Find_Column_to_Ftemove(IOUT), which returns the intrinsic index 
IOUT (associated with the extrinsic index NX, say) of a column located in 
region (i) which may be removed from Fn while allowing the remaining struc- 
ture to be nonsingular. We exchange EPCI, currently located in position IPRI 
of region (Hi), with NX. currently located in position IOUT of region (?). This 
restores the row structure of region (Hi), but EPCI is misplaced in region (i). 
Therefore, we exchange EPCI in position IOUT of region (i) with the extrin- 
sic index located in position MI\C of region (i). By then decrementing MKC, 
we restore the structure of region (?) by decreasing its dimension and that 
of region (?’?) by increasing its dimension. The exchange of NX in I IN with 
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EPCI in IPRI implies a row exchange in To perform this exchange in 

the representation of A^ , we must replace row EPRI of A „ with the current 
representation of row NX of the tableau. Since this current representation is 
not available, we invoke Generate_Row(IIN) to compute it. We then update 
our representation of A j^ 1 by invoking Replace_Row(EPRI,NX). 

9. Pivotal coordinate ((n't), (jjj)). EPCI, the extrinsic index of a basic nonneg- 
ativity constraint currently located in position IPCI of E, is exchanged with 
EPRI, the extrinsic index of a nonbasic nonnegativity constraint currently lo- 
cated in position IPRI of D. The primary exchange properly places EPRI and 
EPCI, so no secondary exchanges are necessary. Since F u is unaffected, no 
additional action is required. 

10. Pivotal coordinate ((?v), (j)). EPCI, the extrinsic index of a basic .E-type 
constraint currently located in position IPCI of B, is exchanged with EPRI, the 
extrinsic index of a nonbasic .E-type constraint currently located in position 
IPRI of D. The primary exchange properly places EPRI and EPCI, so no 
secondary exchanges are necessary. Eu is unaffected, so no additional action 
is required. 

11. Pivotal coordinate ((fu), (jj))- EPCI, the extrinsic index of a basic F-type 
constraint currently located in position IPCI of B , is exchanged with EPRI, 
the extrinsic index of a nonbasic E-type constraint currently located in po- 
sition IPRI of D. The primary exchange misplaces both EPRI and EPCI, 
so both secondary row and column exchanges are necessary. EPRI, cur- 
rently in position IPCI of region (jj), is exchanged with the extrinsic in- 
dex located at position (NXR+1) of region (jj). The column structure is 
restored by incrementing NXR, increasing the dimension of region (;) while 
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decreasing that of region (jj). Note that since the dimension of region (j) 
increases, the dimension of will increase also. Also, the removal of EPCI from 
F n implies that the dimension of F n will decrease. Thus, a column of F n 
must be identified for removal as well. To find such a column, we invoke 
Find_Column_to_Remove(IOUT), which returns the intrinsic index IOUT 
(associated with the extrinsic index NX, say) of a column whose removal from 
F n allows the remaining structure of F n to be nonsingular. We perform a ter- 
tiary exchange by exchanging NX, currently located in position IOUT of region 
(i), with the extrinsic index currently located in position MKC of region (i). 
The structure of region (i) is then restored by decrementing MKC, reducing 
the dimension of region (i) while increasing that of region (ii). To properly po- 
sition EPCI and restore the structure of region (it), we exchange NX, currently 
located in position (MKC+1) of region (ii), with EPCI, currently located in 
position IPRI of region (in). This leaves only NX misplaced. We thus exchange 
NX, currently located in position IPRI of region (in), with the extrinsic index 
located in position MXR of region (ii). The row update is completed by incre- 
menting MXR, increasing the dimension of region (Hi) while decreasing that 
of region (in). Since the dimension of Aj'j 1 has been increased, we must up- 
date our representation. We are required to add the current representation of 
column EPR1, which is available as the pivot column in ROWCOLQ. We also 
require the current representation of row NX of the tableau, which is not cur- 
rently available. We thus invoke Generate_Row(IOUT) With these two rep- 
resentations. we may then invoke Add.Row(NX) and Add_Column(EPCI) 
to complete the update of the representation of Aj", 1 . 
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12. Pivotal coordinate ((iv),(jjj)). EPCI, the extrinsic index of a basic nonneg- 
ativity constraint currently located in position IPCI of B, is exchanged with 
EPRI, the extrinsic index of a nonbasic .E-type constraint currently located in 
position IPRI of D. The primary exchange misplaces both EPRI and EPCI. 
To restore the column factorization, we first exchange EPRI, located in posi- 
tion IPCI of region (jjj), with the extrinsic index located in position (NFR+1) 
of region (jjj). We then exchange EPCI, now in position (NFR-l-1) of region 
(jjj ) , with the extrinsic index located in position (NXR+1). We complete 
the column update by incrementing both NXR and NFR, which has the effect 
of increasing the dimension of region (j), shifting region (jj ) one position to 
the right and decreasing the dimension of region (jjj)- Since the dimension of 
region (j) has been increased, the dimension of Ay will increase as well. To 
restore the row factorization, we exchange EPCI, currently in position IPRI 
of region (iu), with the extrinsic index located in position (MXR+1) of re- 
gion ( iv ). The row update is completed by incrementing MXR, which has 
the effect of increasing the dimension of region (in) while decreasing that of 
region ( 217 ). The dimension of Af, 1 has been increased through the addition of 
column EPCI and row EPRI. Since IPCI is the pivot column and IPRI is the 
pivot row, each is already available in ROWCOLQ. Thus we simply invoke 
Add_Row(EPRI) followed by Add_Column(EPCI) to complete our update 
of the representation of Afj 1 . 
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TABLE 5.2: Secondary and Tertiary Tableau Exchanges 
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VI. FACTORIZATION OF GENERALIZED 
UPPER BOUND ROWS 



A. INTRODUCTION 

We now develop the first of three specializations of the factorization approach. 
We are still interested in the problem: 

(GUB) min : wy 

s.t. : Fy < b F p by n 

Ey < r E m by n 
—Iy <0—7 n by n 

where F, F, —7, w, b and r are as before. We now require that the F-type constraints 
are generalized upper bound (GUB) constraints. Define S’,-, i = 1, . . . ,p to be pair- 

P 

wise disjoint subsets of the set N = {l,...,n} and further define S 0 = N — (J S,- 

i'=i 

( S 0 may be empty). Then S.-fl^t' = 0 for 0 < i < p, 0 < i' < p , i ^ i' and 

p 

(J Si = N. GUB constraints are of the form 
1=0 

* = 1, • • • ,P, hj = ±1 (6.1) 

jes, 

The sets 5,-, i = 0 are called GUB sets and i is the GUB set index. S 0 
identifies those variables that have no coefficient in any GUB constraint. The E- 
type constraints are of arbitrary form. 

A specialization of the Simplex Method to handle GUB constraints was first 
developed by Dantzig and Van Slyke [1967]. They introduce a specialization of 
(GUB) with 8,j = 1 and strict equality constraints. Their algorithm requires a 
working basis of dimension (m + 1). McBride [1972] develops a specialization of 
the mutual primal-dual method to solve a second variant of (GUB) with <5, ; = ±1 
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and strict equality constraints. A dynamic working basis whose dimension cannot 
exceed m is required. A computational analysis by McBride [1972] predicts that the 
performance of this algorithm should exceed that of the Simplex Method when the 
proportion of GUB constraints in the model exceeds approximately 29%, and his 
empirical results tend to confirm this analysis. 

We extend this work by developing a specialization of the factorization ap- 
proach which does not require all GUB constraints to be binding. We will see that 
there are computational efficiencies to be gained by this approach. We first present 
the factored tableau for (GUB). We then discuss the important computational is- 
sues. 



B. THE FACTORED TABLEAU 

Recall from Chapter 4 the primal row basis at an arbitrary point in the solution 
process is: 



B - 





k 


i 


n-l-k 








\ 


O') 


E n 


F 12 


F13 


in) 


F u 


f 12 


f 13 


(jjj) 


V o 


0 


-I J 



}/ 

}k 

} n — l — k 



(6.2) 



We saw in Chapter 4 that it is always possible to identify among the columns of 
B a set of k columns such that the submatrix F n is nonsingular. Since the rows 
of [F n F 12 F 13 ] correspond to GUB constraints, their form is similar to (6.1), with 
column permutations accounting for the difference. Each column of [F n F u Fi 3 ] is 
either zero, or a signed unit vector, and thus by row permutations of these F-type 
constraints, F\\ may be placed into the form of a k by k signed identity matrix An. 
The resulting (GUB) row basis is: 
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CO 


. 

Eu 


E 12 


\ 

E 13 


}l 




(jj) 


An 


Fn 


F, 3 


} k 


(6.3) 


(jjj) 


o 


0 


-I ) 


} n — l — k 





B = 



The corresponding nonbasic constraint rows are then: 



n-l — k 



D = 



(0 


-I 


0 


0 


) k 


(u) 


0 


F 22 


F 23 


}p-k 


(iii) 


0 


-7 


0 


}l 


(iv) 


\ f 2 i 


F 22 


■£-23 j 


} m — l 



(6.4) 



Recall that we wish to observe a one-to-one association between the variables 
whose nonnegativity constraints (which appear in region (i) of D ) are nonbasic 
(nonbinding) and the F-type constraints in region (jj) of B. This association is 
exactly analogous to the idea of “key” variables proposed by Dantzig and Van Slvke 
[1967]. We call the variables whose nonnegativity constraints appear in region (i) 
of D “key” variables, and there is one such variable for each currently basic F- 
type constraint appearing in region (jj) of B. The variables whose nonnegativity 
constraints appear in region (iii) of D are then called “non- key”. 

Using (6.3) and (6.4) and recalling the expression for the explicit transforma- 
tion kernel: 



A„' = (E u - EnF-'Fn )- 1 = (E n - E u ± u F n )-' 
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the principal part of the factored tableau is: 



0) 



(0 

(») 

(iii) 



— An/ii^u' 



~{En — £*i AuFnMn 1 



Uj) 

An 4* AnFij/lii E it An 
FnAl i EiiAn 
— /A 7 i l £n*^u 
-£iiAu + (En — 



Au^i-S — All^l lJ/^i/(£l3 — £l,Ai,F l3 ) 

Fn - " ^iiAuFis) 

i47i'(£i3 *" -CiiAuFis) 

Ejz — £n A h£i 5 — £w/in l (£i3 ” £iiAnFu) 
+£ji AuAj^u 1 (£13 — £u AuFis) / 



(6.5) 

C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 

To examine the actions required by Generate_Column(LC), suppose we 
want to generate column LC of the principal part of the tableau (6.5) and place 
the results in the vector z T = ( 2 j, z 2 , z 3 , z 4 ) T which is partitioned to conform to 
(6.5). Rewriting (6.5) in a more convenient manner, we have: 



9S 



U) 



0 i) 



(0 

(») 



— AiiTu |j4 n ‘ j 

—F T 2 [Arf] 



(i»») 

(iv) —£n [An] ” £ji )J 



-An/*,, [-^u An] + An 

—Fn {”^n ^ii An] 

— Afi^ll An 
-£xi [-^n^i.Aii] 

— £ii An/* i j( — An) + An] 



Uj>) 

— All /'ll |^ll l (^15 £ll Aii/|j)] + An ^13 
-/*i [Aufru - rnAiiFw)] + 



A7i‘(£i3- E\l Alibis) 

— En [> 4 f| l (£i> ~ AiAii/'u)] 
[-AnfutAr/tru- AiA,ir u )) + A„r w ] + £») 



( 6 . 6 ) 

To highlight those computations involving terms consisting entirely of zeros and plus 
and minus ones, vve introduce the following notation: a general matrix product is 
denoted by as in ' ■C-13- A “simple” matrix product (i.e., one in which one 
of the terms consists entirely of zeros, plus ones and minus ones) is denoted by o, as 
in E u o F u - The calculations may then proceed as follows (note that e LC is a unit 
vector with a plus one in position LC): 

1 . Calculate region (Hi): 

(a) If column LC is in (;), then: 

z 3 = A; i i oe LC = (A^) LC 

(b) If column LC is in (jj), then: 

~3 = — An' • (£n) ° (An) /r 



99 



(c) If column LC is in (jjj), then: 



23 = in 1 • \(Eiz) LC - En o An o (F 13 ) LC J 

Notice that {F\ 3 ) LC is either null or it is a signed unit vector, say in row 
k, in which case E n o A n o (F 13 ) LC = ±(E n ) k . 

2. Calculate region (i): 

(a) If column LC is in ( j ), then: 



Z\ — —An o Fi2 o 2 3 



(b) If column LC is in {jj), then: 



z \ — — An 0 Fj2 o 2 3 + (An) tc 

(c) If column LC is in (jjj). then: 



z\ — — An o F12 o z 3 + An o (Ei 3 ) lc . 



Notice that the computation — Fi 2 0 2 3 involves only additions and sub- 
tractions, and that each case in region (i) differs only by an addition 
suffix. 

3. Calculate region (?’?): 
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(a) If column LC is in ( j ), then: 



z 2 — F22 o z 3 



(b) If column LC is in (jj), then: 



z 2 — —F22 o z 3 



(c) If column LC is in (jjj), then: 



Z2 — —F22 O Z3 -{■ (F2z) LC 



4 . Calculate region (tu): 



(a) If column LC is in ( j ), then: 



24 — — -£-21 - z \ ~ E22 • z 3 



(b) If column LC is in (jj). then: 



z 4 



— E 2 1 • Z\ — E22 ’ z 3 



(c) If column LC is in (jjj). then: 



z 4 = -E. 



21 



— E22 ‘ z 3 + (E23) 



LC 
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Then column LC of (6.6) is available a s z T = (zj, z 2, Z3, z 4 ) r . Notice that 
floating point multiplications and/or divisions are necessary only for computation of 
regions (in) and (iv); also note that regions (i), (it) and ( iv ) are linear combinations 
of region (in) and one another, and that each term in these regions differs only by 
an additive suffix. 

D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 

Now we examine the actions of Generate_Row(LR). Suppose we desire to 
generate row LR and place the results of the computation in the vector z = (z 5 ,£ 6 , £7) 
which is partitioned to conform to (6.5). Rewriting (6.5) in a convenient form: 





0 ) 

( ^ . 

n- 4 jy 


0;1 


(jjj) 

J E ,3 [- { j“A n /' ijAn 1 E11 + A n j £13 


(«•) 


— { j”Ai 1 £ ] n^ii'] £n ~ A u 


(»■) 




- { \~F : £11 ] An 


[-^ 32 ^u] ^13 + [- J £n} A n ] £12 Fn 


(««i) 


Mi 


— { £n } An 


M £n ^ \ ~ { l^u 1 ) Em}±u}F u 


(«*) 




- [ [( £?i An£i: ~ AiMn 1 ] £n + £21} An 


[(EjjAjjF ij - EjjMIY! £„ 

* [“ {((£jjAn£l2 - E 3 3)/i7i l ] £ll + £71} All] £53 ~ £53 j 



(6.7) 
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Then we may proceed as follows (note that tm is a unit vector with a plus one in 
position LR): 

1 . Calculate region (j): 

(a) If row LR is in (t), then: 



Z5 — [— (An)£,fl o F12] o A n x 

(b) If row LR is in (it), then: 



?S = (~F22)lR O ^ ii 1 

(c) If row LR is in (Hi), then: 

h = €lr o = (A^)lr 

(d) If row LR is in (iv), then: 

h = {(E 2 i)lr o An o F n - (E 2 2 )lr] • ^n 1 

Notice that the computation of the term (E 2 \)lr o An o — (E 2 2)lr 
involves only additions and subtractions. 

2 . Calculate region (jj): 

(a) If row LR is in (?), then: 

Z(, = (~Zh • E\\ + £-ln) o An 
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(b) If row LR is in (ii), then: 



zq — (—25 • E \ i ) ° A11 



(c) If row LR is in (Hi), then: 



26 — ( — 25 • E n ) O An 



(d) If row LR is in (iv), then: 



26 = (— 25 • Eu — (E2 i)lr) An 



3 . Calculate region (jjj): 



(a) If row LR is in (?), then: 



27 — 25 • E13 + 26 O F 13 



(b) If row LR is in (ii), then: 



27 = 25 • E13 -f Ze O F13 + (F2z)lR 



(c) If row LR is in (212), then: 



27 — 25 • Eiz + 26 o Fiz 
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(d) If row LR is in (tv), then: 



Z 7 = Z S ■ Eiz + 2 6 O F . 13 + (E 22 )lR- 

The row LR of (6.7) is available as z = ( 25 , 26 , 27 ). Note that floating point 
multiplications and/or divisions are avoided for terms involving F n , F 2 2 and F 2 3 , 
that regions ( jj ) and (jjj) are hnear combinations of region (/), and that many 
cases differ only by additive infix terms. 

This computational scheme extends the approach of the mutual primal-dual 
method by introducing additional linear dependencies into the tableau and exploits 
those dependencies to improve efficiency. Empirical evidence suggests that this 
approach also improves numerical stability (McBride [1989]). 

E. DATA STRUCTURES 

The essential information contained in the factored kernel when the factored 
rows are GUB constraints is simply the unique mapping between each basic factored 
constraint and its corresponding “key” nonbasic variable. We require a represen- 
tation of this mapping which is compact, can be efficiently accessed and is easily 
updated. 

One possibility is to maintain the intrinsic ordering of the tableau in such a way 
that the factored row/key variable relationship can be derived implicitly. Region 
(jj) of the tableau corresponds to the basic factored constraints of B. Region (t) 
of the tableau corresponds to the nonnegativity constraints associated with the key 
variables. If we maintain the ordering of (jj) and (?) so that the k th position of (?) 
contains the index of the key variable associated with the factored constraint whose 
index is in the k th position of (jj), we have an implicit representation of Fn. 
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Two difficulties arise with such an approach which lead us to choose an alter- 
nate representation. First, such an ordering complicates the tableau update scheme 
presented in Chapter 5 and potentially increases the computational expense of all 
tableau updates which affect Fn- Second, several update cases require the iden- 
tification of GUB set membership for variables which are not currently assigned 
to region (i) and thus are not currently part of the Fn structure. This may be 
done by accessing the original problem data by column and scanning that column 
for the index of the GUB row, if any, in which the column has a nonzero element. 
Since the original problem data is stored in a super-sparse form, this scheme re- 
quires indirect computer memory addressing. We prefer a method which requires 
less computational expense. 

We thus introduce an additional data structure to manage the structure of 
Fn- KEY(IJ) is an array of type INTEGER, dimensioned from 1 to (m + n), 
which for each basic factored constraint records the extrinsic index of the associated 
key variable, and for each nonbasic key variable records the extrinsic index of the 
corresponding basic factored constraint. Additionally, for each variable which is not 
key (i.e., either non-key variables or basic variables), KEY(IJ) records the extrinsic 
index of the factored constraint for which it may serve as a key variable. If we 
interpret the extrinsic index of the factored constraints as GUB set indices, then 
for every variable KEY(IJ) contains the GUB set index. Note that since GUB 
set membership is fixed for a given variable, all column indices KEY(IJ) can be 
initialized once and need never be updated. Only the factored constraint indices 
require updating. 

F. FACTORED KERNEL UPDATE ACTIONS 

Recall from Chapter 5 the details of the factored kernel update actions 
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Is_Factored_Kernel_Singular, 
Find_Column_To_Remove(/Fii'C), 
Find_Column_To_Add(/FFC),and 
Update JFactored JKernel, 

are factorization-specific. We now discuss these details for the GUB factorization. 

The row-permuted diagonal structure of Fu (An) implies that the question of 
singularity arising from rank-one updates can be resolved strictly on the basis of local 
information. That is, the exchange, deletion or addition of a row leads immediately 
to the identification of a unique GUB set. The search for a corresponding column 
(variable) to complete the exchange, deletion or addition is then limited to those 
sharing this GUB set membership. The simplicity of this structure greatly enhances 
the efficiency of the necessary update actions. 

For example, consider pivot update case ((n), (jj)), in which basic factored 
constraint EPCI, located in tableau position IPCI is removed from the row basis 
and nonbasic factored constraint EPRI, located in tableau position IPRI replaces it. 
Because EPCI was basic, there is a key variable, say EPCIKV, located in position 
IPCIKV of region (i) of the tableau. Because GUB sets are, by definition, pair- 
wise disjoint, EPCIKV may not serve as a key variable for the new basic factored 
constraint EPRI. Thus, a replacement variable among those in region (in) of the 
tableau must be found. Such a replacement must exist, for otherwise B is singular. 
To locate such a variable, we scan the variables in region (in), searching for one 
which is a member of GUB set EPRI. The test is simply this: for each variable JC 
in region (in), if 

KEY (JC) = EPRI. 
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then accept JC as the key variable for EPRI and proceed with the secondary tableau 
update exchange. Otherwise, continue the scan. 

Thus, the action Find_Column_To_A.dd(EKV) scans the indices of region 
(m), performing the GUB set membership test on each index. The first such index 
for which the test is true is chosen as EKV. 

The action Find_Column_To_Remove(EKV) is required when a constraint 
ERI is removed from F n . EKV is then the key variable corresponding to ERI and 
is determined by: 



EKV = KEY(ERI). 

Is_Factored_Kernel_Singular is required when the constraint EPRI is to be 
added to F n and the pivot column index EPCI corresponds to a variable. Thus, 
EPCI is a possible key variable for EPRI and is a convenient candidate. It may be 
accepted as the key variable for EPRI if and only if Fn remains nonsingular after the 
addition of both EPRI and EPCI, and this is the case if EPCI is a member of GUB 
set EPRI. Thus, Is_Factored_Kernel_Singular simply compares KEY(EPCI) and 
EPRI. If they are unequal, F u is singular and the result “true” is returned. Other- 
wise, the result is “false”. 

Finally, we consider Update_Factored_Kernel. Since the GUB set mem- 
bership of a variable does not change, only constraint index updates of KEY(IR) 
are required. For example, if the old key variable index for constraint EPRI is 
EPRIKVO and the new index is EPRIKVN, the required update is: 

KEY(EPRI) = EPRIKVN. 
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VII. FACTORIZATION OF PURE NETWORK 

ROWS 



We remain interested in the problem: 

(PNSC) min : wy 

s.t. : Fy < 6 ; F p by n 

Ey < r ; E m by m 
—Iy <0 ; —I n by m 

where F, E , —I, w, b and r are as before. We now require that each column 
of F consist of zero, one or two nonzero elements. If a column contains a single 
nonzero element within the rows of F, it may be either a plus one or a minus 
one. If a column contains two nonzero elements within the rows of F, one must 
be a plus one and the other a minus one. F-type constraints may also be in- 
terpreted as defining the node-arc incidence matrix of a directed graph. Define 
J = {0, . . . ,p} ,J = {1, . . . ,n} , and A r = {n,-,i G J — {0}} to be a set of nodes and 
A = {a*, k G J | a k = (n,-, ny), i G l,j € j} to be a set of arcs with ordered pairs 
of nodes (tail, head) as elements indexed by k. Note that we interpret node 0 as a 
null node, and thus an arc incident to node 0 is viewed as a single-ended arc. Then 
a graph is defined a s Q = {A r , A). The node-arc incidence matrix of Q is a matrix 
F consisting of p rows (one for each node in A r ) and n columns (one for each arc in 
*4) with elements: 



A. INTRODUCTION 
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Thus, if arc a k = (n,-,7ij) is represented by column k of F , then 



a = 



i th row 



row 



0 \ 

1 
0 

-1 

0 

oy 



and each column of F is either all zeros, contains exactly one nonzero element (which 
may be either plus one or minus one) or contains exactly two nonzero elements (a 
plus one and a minus one). 

A specialization of (PNSC) that has been widely studied arises when the F- 
type constraints are taken to be equalities and the F-type constraints are vacuous, 
resulting in: 



(PNE) nun : 
s.t. : 



wy 

Fy = b 

-Iy < 0 



F p by n 
— I n by n 



Very efficient specializations of the primal Simplex Method have been developed 
and implemented to solve (PNE) (e.g., Srinivasan and Thompson [1973], Glover, et. 
al. [(1974], Bradley, Brown and Graves [1977]). These algorithms exploit some well 
known properties of F (we assume that one redundant row has been removed from 
F): 

1. Every primal Simplex basis B of F (consisting of (p — 1) linearly independent 
columns of F) can be triangulated by row and column permutations, 
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2. Every such basis B is itself the node-arc incidence matrix of a subgraph of Q, 
and this subgraph is a rooted spanning tree, and 

3. F is totally unimodular, implying that if 61 and b? are integer ( p — 1)- vectors 
and xi and x 2 are (p — l)-vectors of unknowns, then for every primal Simplex 
basis B of F the solutions of Bx\ = b\ and x^B = b? are also integer. 

Property (1) allows very efficient execution of Simplex iterations, property (2) mo- 
tivates the use of special data structures which allow efficient storage and update of 
Simplex bases and property (3) allows all computations to be performed in integer 
arithmetic (assuming w and b are integer). 

Several contributions have been made to solving variations of (PNSC) when 
the E-type constraints are not vacuous. These approaches have their roots in work 
by Kaul [1965] and Bennett [1966]. Chen and Saigal considered a version of (PNSC) 
with all equality constraints. Hartman and Lasdon [1972] considered a specialization 
of (PNSC) to the multicommodity capacitated transshipment problem (MCTP), in 
which the E-type constraints are generalized upper bound (GUB) constraints and 
the E-type constraints decouple into independent pure network subproblems. In 
each of these treatments, the resulting algorithm is a specialization of the primal 
Simplex Method. McBride [1972] considered (PNSC) in the context of the mutual 
primal-dual method, with all constraints assumed to be equalities. Each of these 
approaches allows a basis representation which may be factored into two compo- 
nents; an explicit part of dimension m by m and a factored part of dimension p by 
?■ 

Our inequality formulation of (PNSC) allows a dynamic basis representation 
where, just as in the GUB specialization, both the dimension of the factored part 
and the dimension of the explicit part may vary from one iteration to the next. 
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B. THE FACTORED TABLEAU 



Recall from Chapter 3 the primal row basis at an arbitrary point in the solution 
process is: 



(i) 

B = {jj) 

(jjj) 

The corresponding nonbasic constraint rows are then: 




D = 



(0 

(n) 

(m) 

{iv) 



-I 
F 21 
0 

\ F21 



0 

F 22 
-I 
F 22 



and the principal part of the factored tableau is: 



0 \ 
F23 
0 

F23 / 



(7.1) 



(7.2) 





(;) 


(;;) 


0 ;;) 








^ \ 


(>) 


-F n'FnAj 


F,-, 1 + F-'FiAu'EnF^' 


Fu'F xi - FjF u Aj(E X i - E U F^F X2 ) 


(“') 


-FjjA~' + F11F n'Fn/in* 


FjjA^' E\iF u ' — F?\F 


F23 — Fi\ F xx F \3 “f F 2 x F u l F X 2 A IX { E \3 — E XX F U F\z) 






— FiiF^ F u Aii E\\F ,|* 


~~F 22 An( £13 ~~ E u F 1/^3) 
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(7.3) 

The F-type constraints represent the node-arc incidence matrix of a graph. For 
notational convenience, let us define the corresponding directed graph explicitly as 
Q = {Af, A) where (?, A r and A are as previously defined. Note that while we choose 
to interpret all columns of (PNSC) as arcs in C?, any number of these columns may be 



null with respect to the F- type constraints and thus may be interpreted graphically 
as null arcs. 

Since the rows of [Fu F12 F13] in the primal row basis B form a subset of the 
rows of F, [F n F12 F13] may itself be interpreted as the node-arc incidence matrix 
of a graph Qb = {Mb-,Ab}- Mb is a subset of M, but in general Ab - 4 . To 
see this, define Md — M — Mb and consider an arc a q = (n,, n } ) where n,- € Mb 
and rij 6 Md- Then a q “spans” the partitioning of M into Mb and Md- While a q 
is a doubleton arc with respect to Q, it is a singleton arc with respect to Qb- We 
call such an arc a “dynamic singleton”, and we will discover that such arcs require 
special handling in our implementation. 

Since a nonsingular Fn must exist, it is well known that the columns of F n 
represent a rooted spanning tree defined over the nodes of Mb, and that F n may be 
placed into upper triangular form by row and column permutations. If we choose to 
represent F n rather than F {[ 1 in our implementation, we are interested in performing 
two fundamental operations: 

1. Solving linear systems of the type 



F u zi = bi 



and 




where Z\ and z-i are unknown and bi and 62 are rationals (not necessarily 
integers), and 



113 



2. Placing F n in upper triangular form, where F n results from a column ex- 
change, a row exchange, a column and row addition or a column and row 
deletion performed on F n . 

The literature on (PNE) demonstrates that with a proper choice of data structures, 
operation (1) may be performed very efficiently when iq and b 2 are integer, and that 
(2) may be done efficiently when the class of updates is limited to column exchanges. 
We will extend these existing approaches to deal with (1) and (2) in the (PNSC) 
setting. 

C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 

To examine the actions required by Generate_Column(LC), suppose we 
want to generate column LC of the principal part of the tableau (7.3) and place 
the results in the vector z T = (zj, z 2 , z 3 , z 4 ) r which is partitioned to conform to 
(7.3). Rewriting (7.3) in a more convenient manner, we have: 



0 ) 0 ;;] 
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-Fn'Fn Ur,’] 


-Fn'Fn [-Au'EuFF, + FTi' 
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— Eh — En [ — 


—En t —Aii EnF i, 1 ] 
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\ 


-En [-F:i'Fii(-A^E n Fi-i')A F^\ 


-En (£» - EnFTi'Fn)) * Fn'Fn] ) 



(7.4) 
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The calculations may then proceed as follows: 



1. Calculate region (Hi): 

(a) If column LC is in (j), then: 

, 3 = AZ le LC = (Aj) LC 

(b) If column LC is in ( jj ), then: 

(c) If column LC is in (jjj), then: 

*> = ^r, 1 [(£.3 ) 1c - £i.*r, l (F») lc ] 

notice that (F^ 1 )^ and fix' ( Fi 3 ) LC may be computed by backpath traversal 
on F\\ (e.g., Bradley, Brown and Graves [1977] , p. 11). 

2. Calculate region (i): 

(a) If column LC is in (j), then: 

* 11*1 = * 12*3 

(b) If column LC is in (jj), then: 

*ii*i = *12*3 + eLC 

(c) If column LC is in (jjj), then: 



* 11*1 = * 12*3 + (Fiz) LC 
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3 . Calculate region (ii): 



(a) If column LC is in (j), then: 

Zi — —F22Z3 — F21Z1 

(b) If column LC is in (jj), then: 

z 2 = -F22Z3 - F 2 \Z\ 



(c) If column LC is in (jjj), then: 

Z2 = —F22Z3 — F 2 \Zi + (F 2 3) LC 



4 . Calculate region ( iv ): 

(a) If column LC is in (j), then: 

Z\ = -E22Z3 ~ E21Z1 

(b) If column LC is in (jj), then: 

z 4 = —E22Z3 — E 2 \Zi 



(c) If column LC is in (jjj), then: 

Z A = —E22Z3 — E21Z1 + (E 2 3 ) LC 

Then column LC of ( 7 . 4 ) is available as z T = (z l5 z 2 , Z3,z 4 ) T . Note that floating 
point multiplications and/or divisions are necessary only for computation of regions 
(Hi) and (iv), and that regions (i), (ii) and (iv) are linear combinations of region 
(Hi) and one another. 
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D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 



EXPRESSION IN ROW GENERATION 

Now we examine the actions of Generate_Row(LR). Suppose we desire to 
generate row LR and place the results of the computation in the vector z = (£ 5 , z$, zj) 
which is partitioned to conform to (7.3). Rewriting (7.3) in a convenient form: 
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(7.5) 



Then we may proceed as follows: 



1 . Calculate region (j): 



(a) If row LR is in (z), then: 



h = -(Fu')LRF„(Aj) 



(b) If row LR is in (zz), then: 



-1 



s = (F 2 \F n F 12 - F 2 i)lrA n 
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(c) If row LR is in (m), then: 



^5 = (^11 )lr 



(d) If row LR is in ( iv ), then: 



— (-^21 ^* 11 ^ E \ 2 — Ei 22 ) lR ^\\ 



2. Calculate region (jj): 



(a) If row LR is in (i), then: 



•^6-^11 — ~hEn ~ £lr 



(b) If row LR is in (?i), then: 



•^ 6^11 — — (E 2 i)lr 



(c) If row LR is in (Hi), then: 



^6^11 — — -^ 5 -^ 1 1 



(d) If row LR is in (iv), then: 



^ 6^11 — —z$Ei\ — (E 2 i)lr 



3. Calculate region (jjj): 



(a) If row LR is in (i), then: 



27 — Z 5 E\3 + Z&F\3 



(b) If row LR is in (u), then: 



27 — Z5E13 + ZqF\3 + (F 22 )ir 
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(c) If row LR is in (Hi), then: 



Zj — Z$Ei3 -f ZqFi3 



(d) If row LR is in (iv), then: 



27 — 2 5 F 13 + 2 6 F 13 -f ( E^lr 

The row LR of (7.5) is available as 2 — ( 25 , 26 ) 27 ). Note that floating point 
multiplications and/or divisions are avoided for terms involving Fn, F22F13 and F23, 
and that regions (jj) and (jjj) are linear combinations of region (j). 

E. DATA STRUCTURES 

Several suites of data structures have been proposed for algorithms designed 
to solve (PNE) (Johnson [1966], Srinivasan and Thompson [1973], Glover, Karney, 
Klingman and Napier [1974], Bradley, Brown and Graves [1977]). Our implementa- 
tion is based on the last. 

The purpose of the data structures is to define a graph that represents Fn. 
Such a graph forms a rooted spanning tree defined over the nodes of Nb, which 
we denote as Qfu- Note that because Fn is nonsignular, F n must contain at least 
one singleton column. Since Fn is maintained in triangulated form, we associate 
with each row IR the column JC in which the element appearing on the diagonal of 
IR occurs. This association is analogous to the GUB row/key variable association 
defined for the GUB algorithm. We represent this association with the array KEY(), 
of type INTEGER and dimensioned from 1 to (m + n), which records for each row 
of F n the extrinsic index of the corresponding key variable, and for each column of 
F n the extrinsic index of the corresponding basic factored row. KEY() is undefined 
for other row and column indices. 
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We require knowledge of the row ordering of Fn to define a triangulation. 
P0() is an INTEGER array, dimensioned from 1 to p, which for each row IR of F n 
records the extrinsic index of the row that follows IR in the current triangulation. 
This ordering is referred to as preorder, and the successor of IR is called its preorder 
successor. P0() is undefined for factored rows not currently in Fn. Note that the 
triangulated column ordering of Fn may be deduced from P0() and KEY(). 

P0() and KEY() capture the triangulated row and column ordering informa- 
tion for a representation of Fu but provide no means for interpreting this repre- 
sentation as a rooted spanning tree. The predecessor function p() is a well known 
method for representing trees and thus rooted spanning trees. To define p(ik), we 
associate with each node ik of Gf u the row index of the offdiagonal element in the 
k th column of Fu, when the rows are ordered to correspond to a triangulation of 
Fu- Graphically, p(i) may be iterated recursively to trace the backpath from node 
i to the distinguished root node of Gf u - 

We represent p() by the INTEGER array P(), dimensioned from 1 to p, which 
for each row IR in F n records the extrinsic index of the predecessor of IR. It is 
convenient to keep a record of the sign of the diagonal element of each row of F n - 
We use the sign bit of P() to do so. Assume node IRP is the predecessor of node 
IR in the current representation of Fn- If the diagonal element in row IR is a plus 
one, then P(IR) = IRP , while if the diagonal element in row IR is a minus one, 
P(IR) = —IRP. In the graph Gf u , this may be interpreted as follows: if the actual 
orientation of arc a k = KEY(IR) is the same as that recorded in Gf u , which is 
(IR, IRP), then P(IR) > 0). If the actual orientation of arc a k is opposite that 
recorded in Gf u , then P(I R) < 0. 

The final data structure, D(), is an INTEGER array dimensioned from 1 to p. 
For each node IR of Gf u ( row °f Fu) D(IR) records the depth of node IR, where 
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depth is defined as the number of nodes encountered on the backpath between 
IR and the root node. D() allows updates to be performed on the data structures 
representing Qf u ' n a single pass through the data structures, and allows us to avoid 
unnecessary operations in the solution of linear systems to be discussed shortly. 

F. SOLVING LINEAR SYSTEMS 

It is clear from the discussion of Generate_Row(LR) and Generate.Col- 
umn(LC) that the solution of the linear systems, 



and 



zlF n = £>[ 



F\\Z 2 = 62 

where z x and z 2 are vectors of unknowns and b\ and b 2 are vectors of rationals, are 
crucial steps. The practical value of our implementation depends to a large extent 
on the efficiency with which these systems may be solved, and thus we will consider 
this problem in some detail. 

First, we examine the data structures used to represent the terms. Three data 
structures are used to represent and the right-hand side terms. WORK() is a 
DOUBLE PRECISION array, dimensioned from 1 to p, which is used to hold the 
real values of 6] and 62 (we sequence operations so that at any given moment we are 
interested in either bi or b 2 but not both, so that the same array may be used for 
both). An INTEGER array WORKMK(), also dimensioned from 1 to p, is used as 
a nonzero mask for WORK(). The convention is: 
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WORKMK(K) = 0 if VVORK(K) = 0.0 
WORKMK(K) = 1 IF WORK(K) ^ 0.0 



Finally, the array WORKNZ0, of type INTEGER and also dimensioned from 1 
to p, records the extrinsic column (in i>i) or row (in 62 ) index of each nonzero in 
WORK(). The use of WORKMK0 and WORKNZ0 allows the efficient manage- 
ment of nonzeros in f>i and 6 2 . 

The solution vectors Z\ and z 2 are represented by three analogous arrays. 
ROWCOL0 (for rowcolumn()) is of type DOUBLE PRECISION and is dimensioned 
from 1 to (m + n). It is used to store the representation of a row and column of the 
principal part of the tableau, and thus portions of it are used to store the solutions 
of Equation (7.1). Associated with ROWCOL() is ROWCOLMK(), an INTEGER 
array of conformable dimension, used as a nonzero mask for ROWCOLQ, and the 
INTEGER array ROWCOLNZ0, also of conformable dimension, used to record the 
indices of nonzeros in ROWCOLQ. 

In contrast to the usual situation in Simplex-based approaches to linear pro- 
gramming in which each successive right-hand side of the problem may be computed 
by means of a simple update to the previous right-hand side, in this setting such 
is not the case. Thus, we are not able to derive the solution to the current form 
of Equation (7.1) by simply updating the previous solution. Instead, we must solve 
each system from scratch. 

Each right-hand side is itself the result of a sequence of preliminary computa- 
tions. As each new nonzero element is generated by these computations, its value 
is placed in WORKQ, WORKMKQ is updated and the index of WORKQ in which 
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the nonzero appears in placed in the next available position in WORKNZQ. For ex- 
ample, suppose the first nonzero generated for WORK() is the value 6.5 in position 
38. Then the arrays appear as follows: 

WORK(38) = 6.5 

WORKMK(38) = 1.0 

WORKNZ(l) = 38.0 

A counter for the number of nonzeros in WORKQ is maintained (an INTEGER 
variable called NNZW), so at the conclusion of the preliminary computations we 
may iterate through the nonzeros on the right-hand side in the order in which they 
were generated by: 



WORK(WORKNZ(K)), K = 1, ..., NNZW 

Since our representation of F\\ is in upper triangular form, the obvious ap- 
proach for solving 



F n z 2 = b 2 (7.6) 

is by back substitution. All nonzero elements in Fu are either plus ones or minus 
ones, so only assignments, additions and subtractions are necessary. Assuming that 
the current dimension of F n is k by k, we proceed by considering each row IR in 
turn, IR = 1. We assign the value in WORK(IR) to ROWCOL(KEY(IR)) 

and update by an addition the value in WORK() in the row corresponding to the 
predecessor of row IR, (i.e., WORK(P(IR))). 

This approach requires knowledge of the ordering of the rows of Fu in the 
order of last to first in triangulated form, which is precisely the reverse of our data 
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Figure 7.1: A Triangulated Fn 



structure PO(). To support this approach in our implementation, we could define 
and maintain an additional data structure E0(), recording the “endorder” of Fu. If 
we choose not to include an additional data structure, we could instead invert P0() 
in situ whenever this reverse ordering is required, inverting it again when POQ is 
next required. With either approach, it is convenient to make P0() and E0(), if 
included, circular lists and to add a distinguished “artificial node”, say IRA, which 
is then used to locate both the first node (row) in preorder and the first node in 
endorder. To illustrate, Figure 7.1 displays a triangulated form of Fn with row and 
column labels. 

Figure 7.2 presents the corresponding graph Qf u • By assuming the existence 
of an artificial node IRA, we in effect change our paradigm of Qf u to that shown in 
Figure 7.3. 

We may interpret the backsolve method as a labeling procedure on Qf u - I n 
this context, each nonzero in the right-hand side is interpreted as a supply or demand 
at the corresponding node in Qfh- The problem of solving Equation (7.6) is then 
interpreted as one of finding a set of feasible flows defined on the arcs of Qfh- 
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Figure 7.2: A Basis Graph Qp u 

The backsolve approach to solving Equation (7.1) has the advantage of sim- 
plicity, and when viewed as a labeling algorithm has the intuitive appeal of “visiting” 
each node only once. The disadvantages are the need for either additional storage 
(EOQ) or additional computation (two inversions of P0()). 

The right-hand sides of these problems are invariably very sparse; the density 
is typically 0.3% or less. We are thus motivated to explore alternate approaches to 
solving Equation (7.6). The solution approach just outlined requires the “visiting” of 
every node in Qf u - Although WORKMK0 allows us to perform a simple INTEGER 
comparison to determine if the current node has nonzero supply or demand, we have 
strong empirical evidence that suggests that other approaches are more efficient. It 
is important to keep in mind that this problem is very deeply nested within the 
algorithm structure, and must be solved tens of thousands of times for a typical 
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Figure 7.3: An Alternate Graph Paradigm Qp n 

linear programming problem. Small changes in efficiency here have large effects on 
the overall algorithm efficiency. 

Since the right-hand sides are sparse, an alternative is to consider solving Equa- 
tion (7.6) iteratively as a sequence of subproblems, each consisting of a right-hand 
side containing a single nonzero element and a current partial solution. Iterating 
through the nonzeros, we introduce the supply or demand at the corresponding node 
in Qfu > an( l then adjust the flows of all arcs appearing on the backpath of that node 
as necessary. 

Analyzing the worst-case performance of these two approaches, we see that the 
first approach is linear in the number of nodes in Qf u (rows in Fn), while the second 
approach is quadratic in the number of nonzeros in the right-hand side. For the 
densities typically encountered in the problems we have studied, the crossover point 
at which the performance of the linear algorithm overtakes that of the quadratic is 
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in the range of 400 to 1000 rows in F\\. However, in the testing we have performed, 
the worst-case analysis is pessimistic, and in practice the performance of the second 
approach always dominates that of the first. 

The second linear system of interest is 

zfF u = bl . (7.7) 

We interpret the problem as one of being given flows on the arcs of Qf u and being 
asked to find the necessary supplies and demands at the nodes. Analogies to each of 
the approaches to solving Equation (7.6) may be found for solving Equation (7.7), 
and our empirical evidence strongly suggests that an approach analogous to the 
second method for solving Equation (7.6) is preferable for solving Equation (7.7). 

Having made the design decisions for solving Equation (7.6) and Equation 
(7.7), it is worthwhile to review our paradigm for Qf u - Our preferred solution 
techniques do not require a partial ordering of all rows of F\\. Rather, we require 
partial ordering information only among each set of coupled rows, or, graphically, 
among each disjoint component of 0f u - Rather than introducing an artificial node, 
we may treat each disjoint component as an independent entity. It is then conve- 
nient to treat singleton arcs (either dynamic or static singletons) as self-loops, and 
construct the predecessor function so that it forms a ring within each component. 
Our graph paradigm of Figure 7.1 then becomes as shown in Figure 7.4. The data 
structure representation of Figure 7.4 is then: 
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Figure 7.4: The Implementation Paradigm for C/p,, 

G. FACTORED KERNEL UPDATE 

The update of a rooted spanning tree representation of a pure network tri- 
angulated basis which results from a column exchange is well understood and has 
been treated thoroughly in the literature. Bradley, Brown and Graves [1977] give 
an excellent account of such an update in an algorithmic setting very similar to the 
one here. Hence, we will not repeat the details. In this setting, however, a column 
exchange is but one of four possible updates required. We must also deal with row 
exchanges, row and column additions and row and column deletions. Our general 
approach will be to reduce these three cases to the column exchange case through 
a sequence of operations which may be thought of as pre- and post-processing. The 
following three cases, along with the column exchange case, comprise the actions 
required for Update JFactored_Kernel. 

The first case we consider is the row and column deletion case, in which the 
dimension of Fu decreases. It is convenient to limit row/column deletions to row/key 
variable pairs. This is because removing a row/key variable pair always preserves 
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nonsingularity of the remaining factored kernel, and because, with our choice of data 
structures, given one member of the pair, it is easy to identify the other member 
through the use of the KEY() array. 

With our paradigm for , the key variable for row IR corresponds to the arc 
which connects node IR to its predecessor. The removal of IR and KEY(IR) creates 
a new disjoint component within Qp u consisting of all nodes and the associated key 
variables on whose backpaths IR and KEY(IR) appear. If IR is a leaf (a node which 
is incident to a single arc), then the new disjoint component is null. Within the new 
disjoint component, D() changes for all nodes, but P() and P0() change only for the 
first node in preorder, say IRFIRST, and the last node in preorder, say IRLAST. 
For the first node in preorder, its predecessor becomes itself (with the sign bit used 
to indicate arc orientation), forming a self-loop. The update is thus: 

P (IRFIRST) = IRFIRST 

Its depth is updated to be zero, and the magnitude of this depth change is recorded 
for future use. For every other node in the new disjoint component, Z)() is reduced by 
the magnitude of the IRFIRST depth change. Finally, PO() of IRLAST is changed 
to the first node in preorder, restoring the local ring structure of P0(). The update 
is: 



PO (IRLAST) = IRFIRST 

These updates can be accomplished in one pass through the data structures of the 
new disjoint component. The structure of Qfu is restored by changing P0() of the 
predecessor of IR to the index of the node which was the successor of IRLAST prior 
to the update. 
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The second update case is a row and column addition. Suppose IRADD is the 
row (node) to be added to F u (C?/r n ). We wish to treat this update as a column 
exchange. To do so, we first incorporate IRADD into the structure and associate 
with it an imaginary (“logical”) arc which forms a self- loop. We then perform a 
column exchange with the new column, say JCADD, entering Qp u and the logical 
arc leaving Qf ix - 

The difficulty in incorporating IRADD into lies with the dynamic single- 
ton arcs. Qfu currently may have many singleton arcs which appear in our paradigm 
as self-loops. If any of these dynamic singletons are actually incident to IRADD in 
G, incorporating IRADD into requires changing their representation. 

To illustrate this point, assume IRADD has the node label “9” and we wish 
to incorporate it into the graph shown in Figure 7.4. Assume also that arcs 103 and 
107 are dynamic singletons which are actually incident to node 9; that is, arc 103 
= (9,4) and arc 107 = (9,8). We first initialize the data structures for node 9 as 
follows: 



D{ 9) = 0 
P{ 9) = 9 
PO( 9) = 9, 

which has the effect of placing node 9 at depth 0 as a disjoint component consisting 
of a single node with a (logical) self-loop, as shown in Figure 7.5. 

We then change the representation of arcs 103 and 107 from dynamic singletons to 
doubletons. To do this, we change P(4) from -4 to -9 and P(8) from -8 to -9. We 
then assert a partial ordering among those (formerly) disjoint components which 
have been merged by the change in the representation of the dynamic singletons, 
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Figure 7.5: Initialization of Node 9 Prior to Incorporation into Qp u 

and retriangulate. This is accomplished by increasing the depth of each node in the 
component by one, changing PO(7) from 4 to 8, and P0(8) from 8 to 9. This update 
requires a single pass through the data structures of each affected component. After 
incorporation, Qf u appears as shown in Figure 7.6. We may now complete the 
update by performing a column exchange. 

To identify the dynamic singletons which are affected by such an update, we 
access the extrinsic problem data structures for row (node) IRADD. For each column 
with a nonzero in row IRADD, we test whether or not that column is currently key 
for some row in F\\. If so, that column is currently represented in Qf u as a dynamic 
singleton. 

The final update case is the row exchange case. Suppose we want to replace 
row IROUT with row IRIN. We treat this in two stages. The first stage is a row 
and column deletion case, in which the node IROUT and the arc KEY(IROUT) are 
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removed from Qf u - To complete the update, we seek a column (arc) which may 
be added to Qp^ along with row (node) IRIN. We require an arc which is either 
a static singleton, incident to node IRIN, or a doubleton, incident to node IRIN 
and a node (other than IROUT) which is currently in Gf u . We first consider arc 
KEY(IROUT). If is satisfies the second condition (obviously, it cannot satisfy the 
first condition), we designate it as the arc to be added. If not, we search among the 
variables (arcs) in region (Hi) of the tableau for such an arc. We know one must 
exist, for otherwise B is singular. We select the first such arc found as the arc to be 
added. We then perform the row and column addition case. 

Recall that the action Is_Factored_Kernel_Singular() is required when we 
have identified an arc (column) for removal from Qp u and we are considering a 
candidate arc to replace it. We want to determine if this exchange preserves the 
nonsingularity of Fu. The fact that Qp u is a rooted spanning tree provides the 
necessary structure to support a simple test of nonsingularity. Assume JCOUT is 
the (known) arc to be removed from Qp n and JCIN is the candidate replacement. 
In our paradigm, F n is nonsingular if and only if each disjoint component of Qp 11 
is connected and contains exactly one cycle, which must be a self loop occurring at 
the component’s root (the root is the distinguished node in the component whose 
depth is zero). The removal of JCOUT creates a new disjoint component which has 
no root self loop. If the addition of JCIN fails to correct this, the proposed exchange 
is singular. 

The specific test is this: identify the nodes incident to JCIN using the original 
problem data structures. Note that there may be zero, one or two such nodes. If 
none of these nodes are currently included in the structure of G f 11 (F\i), then the 
proposed exchange is singular. For each such node which is currently included in 
the structure of Gf u , traverse the backpath of that node by recursively iterating 
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the predecessor array P() until either the root node of that component is found, or 
until the “join” is located (the “join” is that node with largest depth which appears 
on the backpath of both nodes incident to JCIN. The join exists only if JCIN is 
incident to two nodes, both nodes are currently included in the structure of Qf U i 
and the two nodes are in the same disjoint component.). If JCIN is encountered 
during this backpath traversal, the proposed exchange is nonsingular. Otherwise, it 
is singular. 

The action Find_Column_to_Remove(IOUT) is particularly simple, since, 
as mentioned previously, we always remove row/key variable pairs. Thus, if IR is 
the node to be removed, then KEY(IR) is the corresponding column to be removed. 

Finally, the action Find_Column_to_Add(IIN) requires searching among 
the variables in region (Hi) of the tableau. For each candidate arc, we invoke 
Is_Factored_Kernel_Singular(). If the response is “FALSE”, we have found the 
column to be added. Otherwise, we continue the search. 
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VIII. FACTORIZATION OF GENERALIZED 

NETWORK ROWS 



A. INTRODUCTION 

The problem of interest remains: 

min : wy 

y 

(GNSC) s.t. : Fy < 6; F p by n 

Ey < r; E m by n 
— Iy < 0; —I n by n 

where F , E , —I, w, b and r are as before. We now require that each column of 
F have at most two nonzero elements, which may be of arbitrary sign. We may 
associate a generalized network with F by defining a node corresponding to each 
row i of F and an arc F J corresponding to each column j of F. The arcs are defined 
as: 




M) 

M) 

( 0 , 0 ) 



if Fij ^ 0, Fkj 0 and i < k < p 
if Fij ^ 0, Fkj = 0 for all k ^ i < p 
if Fij = 0 for all i < p 



We let A — {-F 1 , . . . , F n } ,A r = {1, . .. ,p} and define the graph Q = {A r ,A}. Arcs 
of the form (t, 0) are singleton arcs, sometimes called root arcs. Arcs of the form 
(0, 0) are null. 

The most widely studied specialization of (GNSC) is obtained when the F- 
tvpe constraints are equalities and the £-type constraints are vacuous, resulting 
in: 
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(GNE) 



min : 
y 


wy 




s.t. : 


Fy < 6; 


F p by n 




Ey < r 


E m by n 




-Iy < 0; 


—I n by n 



Dantzig [1963] provides the seminal treatment of (GNE), identifying the im- 
portant structure that leads directly to efficient primal Simplex-based algorithms 
for its solution. Implementations have been reported by Glover, Klingman, Hultz, 
Karney and Elam [1972, 1973, 1977, 1978, 1979] and Brown and McBride [1984]. 

(GNE) may be viewed as a generalization of (PNE), and the cost of such 
generalization is a weakening of the properties which so strongly characterize (PNE) 
and which lead to the efficient and elegant implementations. In (GNE), F is not 
totally unimodular, and thus an implementation may not be restricted to integer 
arithmetic. It is not possible to triangulate every primal Simplex basis extracted 
from F by row and column permutations. The subgraph generated from the columns 
of a primal Simplex basis no longer form a rooted spanning tree defined on the nodes 
of Q. Apparently, much has been given up in the generalization. 

In fact, significant structure remains in (GNE). A well known result, due to 
Dantzig [1963], shows that any primal Simplex basis extracted from F can be put in 
the form of Figure (8.1) by row and column permutations. Each square submatrix 
component B k is either upper triangular or nearly upper triangular with only one 
element below the diagonal. Notice that if each B k is strictly upper triangular, the 
structure is analogous to that found in the (PNE) primal Simplex basis. 

Interpreting the structure of B k as a subgraph of Q, we find that when B k 
is upper triangular, the subgraph may be viewed as a disjoint component with a 
single self- loop at the root node, exactly as in the (PNE) case. When B k is not 
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' B l 
B 2 



B k 



Figure 8.1: Near-Triangulated Simplex Basis Corresponding to (GNE) 

upper triangular, the subgraph still forms a disjoint component with a single loop, 
but that loop is no longer a self-loop. For example, Figure (8.2) shows a nearly 
triangular B k with row and column labels. 
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105 
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Figure 8.2: A Sample Nearly- Triangulated (GNE) Simplex Basis Com- 
ponent 

The resulting subgraph is shown in Figure (8.3). Such a subgraph is commonly 
called a “one-tree”, (an apparent oxymoron). 

This structure allows the extension of the algorithm and data structures de- 
veloped for (PNE), leading to efficient implementations to solve (GNE) (e.g., Brown 
and McBride [1984]). 

(GNSC) has received less attention in the literature than (PNSC). Hultz 
and Klingman [1976] develop a primal Simplex-based approach to an equality- 
constrained formulation of (GNSC), and report on an implementation which allows 
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Figure 8.3: “One-tree” Subgraph 

a single F-type constraint (Hultz and Klingman (1978)). McBride [1985] develops 
an algorithm for solving a generalization of (GNSC) in which complicating columns 
as well as complicating constraints are permitted, and reports on an implementation 
of that algorithm. 

B. THE FACTORED TABLEAU 

The algebraic development of the primal row basis D, the nonbasic rows D 
and the factored tableau is exactly as shown in Chapter 7 and is not repeated here. 
Note that Fu, the factored kernel, may now be placed in the form shown in Figure 
(8.1) by row and column permutations. The corresponding graph, consists of 
one or more disjoint components, each of which contains exactly one loop. The loop 
may be either a self-loop, as in the case of the disjoint components of (PNSC), or 
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it may be a “root loop” (a loop consisting of two or more nodes, all of which are at 
depth zero in the component), as shown in Figure (8.3). 

C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 

The sequencing of computations for column generation in (GNSC) is exactly 
the same as that for (PNSC). 

D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 

The sequencing of computations for row generation in (GNSC) is exactly the 
same as that for (PNSC). 

E. DATA STRUCTURES 

Several proposed data structures for algorithms tailored to solve (GNE) have 
been offered (Glover, Klingman and Stutz [1973, 1974], Elam, Glover and Klingman 
[1979], Brown and McBride [1984]). Our implementation is based on the last, which 
extends to (GNE) the data structures developed in Bradley, Browm and Graves 
[1977] to solve (PNE). 

We have seen the similarity between the structure of the disjoint components 
that arise in the graph corresponding to F n in (GNSC) and the structure of the 
components in (PNSC). To develop a representation of these components, we wish to 
extend the data structures developed for (PNSC) in a natural way. Knowledge of the 
row ordering of Fn is again maintained in an INTEGER array PO(), dimensioned 
from 1 to p. As before, it is convenient to make PO() a circular list defined on 
each disjoint component. The unique correspondence between each row of F n and 
a distinguished column is maintained in the INTEGER array KEY(), dimensioned 
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from 1 to (m + n). The conventions for KEY() are exactly as in (PNSC). The depth 
array, D(), is defined exactly as before. However, as suggested by our diagram 
in Figure ( 8 . 3 ), we adopt the convention of placing all nodes appearing on a loop 
at depth zero, which is consistent with our definition of “root loop”. Finally, a 
representation of the predecessor function is needed, and we define P() as before. 
For any row IR in the factored kernel, P(IR ) > 0 indicates that the diagonal element 
in the near-triangulation is the first non-zero factored element in its column. 

F. SOLVING LINEAR SYSTEMS 

Just as in (PNSC), the crucial steps in the generation of rows and columns of 
the principal part of the tableau are solving the systems 

zjF u =bJ ( 8 . 1 ) 

and 

F \\ z 2 — ^2 (8-2) 

where Z\ and zi are vectors of unknowns and b x and 62 are vectors of rationals. Our 
approach for solving these systems closely parallels that developed for (PNSC). 

The data structures used to represent ^1 , ^2 > ^1 an ^ ^2 are exactly the same 
as in the (PNSC) implementation. To represent b x and b 2 , we use: 

WORK() 

WORKMKQ 

WORKNZQ, 
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and for Z\ and 2 2 ' 

ROWCOL() 

ROWCOLMK() 

ROVVCOLNZ(). 

The right-hand side sparsity-exploiting approach for solving (8.1) and (8.2) 
developed in Chapter 7 may still be used, but two complications arise in the (GNSC) 
setting. First, the nonzero elements of F u are no longer restricted to be plus and 
minus ones, and may assume arbitrary values. Thus, we are not able to restrict all 
arithmetic operations to additions and subtractions. 

Interpreting (8.1) and (8.2) as problems of finding balancing supply and de- 
mands (8.1) or feasible flows (8.2), we see that in (GNSC) we must allow for gains 
and losses in the flows over the arcs of Qp u . Because our formulation allows two arbi- 
trary nonzeros in each column of F (rather than one arbitrary nonzero and one plus 
one, for example), backward and forward substitution schemes require both multi- 
plications and divisions. We may eliminate the need for divisions by pre-computing 
both ratios of the nonzero elements in each arc. That is, if a column has nonzero 
elements a and b in rows of F, we pre-compute the ratios | and £. We then substi- 
tute the pre-computed value for the division whenever it occurs. Storing a pair of 
DOUBLE PRECISION real numbers for every column of (GNSC) is wasteful, and 
instead we define a single pointer for each column which points to the location of 
the first of two DOUBLE PRECISION real values. The first is the value | and the 
second is £. (Note that for singleton arcs, the values a and £ are stored instead.) 
We compute and store only the unique ratio pairs for a given problem instance, and 



140 



the number of such unique pairs is usually quite small: typically, a problem with 
10,000 to 20,000 columns has only 30 to 40 unique ratio pairs. 

The second difficulty is that the disjoint components of Fn may be nearly tri- 
angular rather than triangular, and thus backward and forward substitution cannot 
be applied directly. However, Dantzig [1963] shows that a variation of backward 
and forward substitution may be used to solve (8.1) and (8.2). This method solves 
the triangulated part of a disjoint component exactly as is done in backward and 
forward substitution. When a root loop is encountered, the method requires two 
calculations for each row 7 or column of F\\ corresponding to a node or arc on the 
loop. The second calculation involves a term Dantzig [1963] called the “loop factor”, 
which is common value for every node (or arc) on the loop. Our implementation 
uses this modified forward and backward solution technique for solving (8.1) and 
(8.2), and we store unique values of loop factors in a DOUBLE PRECISION array 
called GNROOT(). 

G. FACTORED KERNEL UPDATE 

We first consider the actions required by Update_Factored_Kernel. 

Brown and McBride [1984] give an excellent treatment of the update required 
for a column exchange within a generalized network basis in an algorithmic setting 
very similar to the one here, and thus we do not repeat that discussion. As in 
the case of (PNSC), however, three additional classes of updates may occur in this 
setting: row exchanges, column and row additions and column and row deletions. 
We again treat these three cases by reducing them to the column exchange case 
through sequences of pre- and post-processing operations. 

We again limit row and column deletions to row/key variable pairs. When 
row IR and its associated key variable KEY(IR) are removed from F n , a disjoint 



141 



component is created which has no loop. If IR and KEY(IR) do not appear on a root 
loop of Q Fni a new disjoint component is created, just as in the (PNSC) algorithm. 
If IR and KEY(IR) do appear on a root loop of Gf u , the resulting retriangulated 
component contains a self-loop rather than a root loop. The update of the subgraph 
data structures is very similar to those presented in Chapter 7. 

Row and column additions require the incorporation of the incoming node into 
the structure of prior to the column exchange, just as in (PNSC). The technique 
for (GNSC) is exactly the same as in (PNSC). Once the new node has been added 
to Qf U i a column exchange operation is performed, with a logical arc (which forms 
a self-loop on the new node being added to Qf u ) being replaced by the new arc 
(column) being added. 

Finally, row exchanges are handled as a two-stage update, exactly as in (PNSC). 

The task of determining the singularity of Fn (Is_Factored_Kernel_Singu- 
lar) is more challenging in (GNSC) than in either of the other two algorithms due 
to the arbitrary nonzero elements. In each of the previous two algorithms, we have 
been able to determine with certainty the singularity of Fn indirectly. In the (GUB) 
factorization, Fn = An, a signed identity matrix. Since An is orthogonal (i.e., 
A?i- An — n, Fn is perfectly conditioned (see e.g., Golub and Van Loan [1983] for 
a discussion of matrix conditioning). In (PNSC), Fn may be triangulated, placing it 
in a form with plus and minus ones on the diagonal. Thus, the determinant of F n is 
either plus or minus one. We may therefore determine its singularity by examining 
the structure of Qf xx rather than considering Fn directly. In either case, as long 
as the accumulated round-off error in the current representation of the problem 
is such that we can distinguish plus or minus one from zero, we may rely on the 
structure of our representation of Qf xx to deduce the singularity of Fn- Thus, we 
are able to discern singularity logically and need not resort to analytic methods. 
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In (GNSC), however, we may not rely exclusively on the structure of Qp n to reach 
conclusions about the singularity of F u . It is not difficult to invent examples of 
factored kernels corresponding to instances of (GNSC) which are ill-conditioned. 
We expect ill-conditioning of Fu to lead to serious numerical difficulties, which 
we wish to avoid if possible. Since we generally have freedom in determining the 
structure of Fn and since our fundamental rank-one update of Fn is the column 
exchange, we use a heuristic based on column exchanges to anticipate conditioning 
problems. Recall in the column exchange update, a column has been selected for 
removal from Fu and a candidate column has been identified as its replacement. 
We use the backpath traversal method described in Chapter 7 to determine if the 
proposed exchange maintains the required structure of Qf u - ^ it does not, we can 
conclude that the exchange renders Fu singular and we reject the candidate arc. 

As we traverse the backpath, we perform calculations which, if the exchange 
does in fact preserve the required structure of Gf u , will ultimately compute the 
determinant of the disjoint component in which the replacement arc will appear if the 
proposed exchange is performed. If the absolute value of the computed determinant 
is too small (a decision controlled by a implementation parameter), we conclude that 
the proposed update produces a submatrix of Fu (the submatrix corresponding to 
the new disjoint component) that may be ill-conditioned, and thus Fn may itself 
be ill-conditioned. We therefore reject the proposed candidate arc. 

This approach is attractive because it is computationally inexpensive and, 
based on our empirical evidence, it works well in practice. It is of course merely 
a heuristic, since it is easily shown (e.g., Golub and Van Loan [1983)] that the 
determinant can be a poor indicator of matrix condition. Further, we do not actually 
compute the determinant of Fu, but rather only the determinant of a submatrix of 
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The remaining update actions, Find_Column_to_Add(ICI) and Find.Col- 
umn_to_Remove(ICl) are treated in a manner analogous to that in (PNSC). 
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IX. COMPUTATIONAL RESULTS 



A. INTRODUCTION 

The algorithms developed here have been implemented and extensively tested 
within the framework of a commercial-quality optimization system: the X-System 
of Brown and Graves [1975]. This system employs the Graves mutual primal-dual 
algorithm in a variety of large scale optimization applications, including linear, non- 
linear, mixed integer and decomposed models. Although we report computational 
results only for linear models, our factorizations are seamlessly integrated into the 
X-System and support all system features. 

B. TEST PROBLEMS 

The benchmark test suite is drawn from a wide variety of actual applications. 
Table (9.1) provides a short synopsis of each model, quoting from the abstract 
where a reference in the open literature is available. In some cases, several different 
instances of models are reported. We selected these models because they provide 
a representative sample of size, structure and taxonomy of contemporary real-life 
applications of linear programming. For those models that are mixed integer, we 
report solution statistics for the initial linear programming relaxation. 

TABLE 9.1: Description of Test Suite Models 

• GTE “The seven Telephone Operating Companies within GTE have adopted 
an integrated business system called Capital Program Management System 
(CPMS) to guide their 3 billion dollar per year capital planning. The system 
includes a large scale mixed integer programming optimization system that 
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optimizes the critical economic tradeoffs between maximizing the long-term 
budget value of the firm’s equity and satisfying shorter-term financial con- 
straints, resource limitations and service objectives. Investment opportunities 
for the next 5 years are modeled as 0-1 variables with alternative implemen- 
tations for each. The objective is to maximize the net present value of the 
capital portfolio. There are financial constraints on capital, internally gen- 
erated funds, net income to common, and limits on resources such as labor 
hours, lines installed, etc. There are also constraints that enforce logical re- 
lationships among opportunities (such as, if choose A then must choose B).” 
See Bradley [1986]. 

• INVEST “Capital allocation and project selection for a large multi-national 
firm is modeled as a two-stage multi-year nonlinear capital budgeting problem 
w'ith over 40,000 integer variables. Innovative modeling yields subproblems 
easy to solve, and optimality is achieved with a single iteration of the nonlinear 
master problem.” See Harrison, Bradley and Brown [1989]. The instance 
reported here is a linear program subproblem of the two-stage model. 

• TANKER “A crude oil tanker scheduling problem faced by a major oil com- 
pany is solved using an elastic set partitioning model. The model takes into 
account all fleet cost components, including the opportunity cost of ship time, 
port and canal charges, demurrage and bunker fuel. The model determines 
optimal speeds for the ships and the best routing of ballast (empty) legs, as 
well as which cargoes to load on controlled ships and which to spot charter. All 
feasible schedules are generated, the cost of each is accurately determined and 
the best set of schedules is selected.” See Brown, Graves and Ronen [1987]. 
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• GAS “Natural gas utilities supply about one fourth of the energy needs of 
the United States. From wellhead to consumer, operations are governed by an 
astounding diversity of purchase, transport and storage contract agreements 
which enable complex physical distribution systems to meet future demands 
no more predictable than next year’s weather. Gas is a highly detailed op- 
timization model which utilities use to plan operations and to justify such 
plans to regulatory agencies is developed.” See Avery, Brown, Rosenkranz 
and Wood [1989]. 

• ODSM “A commonly occurring problem in distribution system design is the 
optimal location of intermediate distribution facilities between plants and cus- 
tomers. A multicommodity capacitated single-period version of this problem 
is formulated as a mixed integer linear program. A solution technique based 
on Benders Decomposition is developed. . . . An essentially optimal solution 
is found and proven with a surprisingly small number of Benders cuts.” See 
Geoffrion and Graves [1974]. The instances reported here are decomposition 
master problems. 

• TAM “The annual decision on how much of the Air Force procurement budget 
should be spent on the many different aircraft and how much should be spent 
on the many different munitions is of great interest to many people. How 
the Air Force staff develops information to support the decision has changed 
over the years. Currently, a linear program is being used by the Air Force 
Center for Studies and Analysis and is being tested by the Munitions Division 
of the Plans and Operations Directorate (AF/XOXFM) for munitions tradeoff 
analysis. The LP uses existing data and estimates on (1) aircraft and munition 
effectiveness, (2) target value, (3) attrition, (4) aircraft and munition costs, 
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and (5) existing inventories of aircraft and munitions. Other factors considered 
are weather and length of the conflict.” See Might [1987] and Jackson [1989]. 

• PHOENIX “The U. S. Army operates a fleet of over 7,000 helicopters to 
perform combat and combat support tasks. Although newer, more technically 
advanced helicopters have been and are being procured, the majority of the 
fleet is still composed of helicopters that were built during the late 1960s 
for use in South Vietnam. These older airframes are rapidly reaching the 
end of their useful lives and must be (i) replaced by newer, more advanced 
designs, (ii) gutted and refit or replaced by a combination of (i) and (ii). Army 
force planners recognize that this problem can only be solved by a long-term 
program and the commitment of billions of dollars. Phoenix is a software 
system employing mixed integer linear programming to help Army aviation 
staff and officers develop long- range helicopter modernization plans.” See 
Clemence, Teufert, Brown and Wood [1988]. 

• AMM04H “A four-commodity transshipment model for delivery over time of 
military products from production and storage locations to overseas locations 
to support theater operations is developed. The model covers five physical 
echelons, including production plants, storage depots, ports of embarkation, 
ports of debarkation and geographic field locations. Road, rail, sea and air 
transportation are modeled, and product demands are time- phased. Capaci- 
tation occurs primarily on sea and air links, and as throughput capacities on 
transfer points, requiring replication of some echelons. The objective of the 
model is to minimize deviation from on-time deliveries.” See Staniec [1984]. 

• GK A weekly multi-plant production/inventory/transshipment linear pro- 
gram from a consumer products industry is developed. The model is meant 
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to guide weekly processing and packaging decisions. Production consists of 
two stages: basic products are produced and then packaged into different- 
sized containers to yield finished products. Processing lines typically produce 
a subset of the basic products and have limited capacity with overtime charges 
for weekend shifts. Packaging Lines for finished products are similar. In-house 
inventory capacity is limited although outside storage is available at additional 
cost. Inter-plant shipments are made by rail or truck. See Wood [1989]. 

C. METHODOLOGY 

We wish to evaluate implementations of our three algorithms on the basis 
of computation time and computer memory requirements. Since each algorithm is 
simplex-based, the formal theory of algorithmic complexity provides no basis for 
preferring one to another, since in the w'orst case none enjoys a measure of running 
time that is polynomial in the size of the problem specification (see e.g., Garey 
and Johnson, [1979] for a discussion of algorithmic complexity and Klee and Minty 
[1972] for an analysis of the Simplex method). Thus, w r e are led to consider “typical” 
performance by gathering empirical evidence on the performance of the algorithms 
solving “typical” problems. 

We prefer an implementation that is both fast and requires little computer 
memory. Most researchers who have reported on implementations of related al- 
gorithms have been concerned primarily with execution speed, and certainly it is 
important. However, we have seen in our algorithmic setting that once high speed 
storage has been allocated for the problem representation and for the program code, 
all remaining memory is available to store the representation of the explicit trans- 
formation kernel and the factored kernel. If the solution trajectory is such that their 
combined size never exceeds available memory, we succeed in solving the problem. 
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If not, we fail. When success or failure depends on the total storage requirements 
of the inverse representation, we may be willing to sacrifice execution speed in ex- 
change for an economical representation of the inverse. This is a classic theme in 
computer science and software engineering, and we believe its importance in this 
context has been largely overlooked. 

We will be comparing the performance of four separate implementations. The 
first is an unadorned version of the X-System, which implements the Graves mutual 
primal-dual method as presented in Chapter 2. There is no underlying factorization. 
We refer to this implementation as “X”. The second implementation is the GUB 
factorization presented in Chapter 6, and referred to as “GUB”. The third is the 
pure network factorization of Chapter 7, referred to as “PNET”, and the last is the 
generalized network factorization, called “GNET”, of Chapter 8. 

The ideal approach for this computational study, would be to develop four 
equivalent formulations of each model, each customized for its particular implemen- 
tation with the goal of inducing a large factored row set of the appropriate type. 
This approach is a consistent theme in the literature dealing with specialized al- 
gorithms and one that we strongly endorse. Alternate formulations of a model are 
often available, and it seems sensible to choose one that as strongly as possible 
exploits the strengths of the solver. 

However, all of the models used here are “off-the-shelf” in the sense that they 
were developed at various times by various modelers, and we cannot afford to develop 
alternate formulations. Thus, our approach is to preserve a single, unfactored rep- 
resentation of each model, and attempt to identify favorable row structures through 
the use of heuristics. Our procedure is based on the work of Brown and Thomen 
[1980], Brown and Wright [1983] and Brown, McBride and Wood [1985]. The heuris- 
tics are greedy and myopic in the sense that they initially consider the entire row 
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set of the problem, and discard one row at a time without backtracking until the 
remaining set satisfies the desired row factorization. This may easily result in the 
confounding or destruction of structure intended or perceived by the modeler. While 
our approach yields interesting and useful observations about the implementations, 
it is in some ways a poor substitute for the customized model method. 

Table (9.2) tabulates the important structural information concerning the 
model instances we will be solving. The column headings may be interpreted as 
follows: m is the total number of structural constraints (note that this use is differ- 
ent from that in the problem specifications of Chapters 6, 7, and 8), n the number 
of variables, Pgub the number of GUB rows identified by the heuristic, pp ^ the 
number of pure network rows, pgn the number of generalized network rows and NZ 
the number of nonzeros in the technological coefficient matrix of the model. For 
example, consider the first model in Table (9.2), GTE. The structural constraints 
contain 57,563 nonzeros, and the model consists of 6,624 variables. When viewed 
as an unfactored mutual primal-dual model, it consists of 960 explicit constraints 
and 0 factored constraints. When viewed as a GUB factorization, it consists of 909 
factored (GUB) constraints and 960 — 909 = 51 explicit constraints. Similarly, when 
viewed as a pure network factorization (PNET), it contains 909 factored rows and 51 
explicit rows. Finally, when viewed as a generalized network factorization (GNET), 
it consists of 922 factored rows and 38 explicit rows. 

D. COMPUTATIONAL RESULTS 

We solved each of these problem instances on an IBM 3033/ AP under the 
VM/CMS operating system using VS FORTRAN 1.4.1. A virtual machine size 
of six megabytes was used, simply because it is the largest size normally available 
to us. It is the nature of a time-shared system that measurements of processing 
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times are somewhat imprecise due to system load factors and accounting techniques 
for system processing overhead. We have attempted to mitigate these effects by 
performing our experiments during periods of low system usage. 

Table (9.3) displays the solution times for each of the test problem instances by 
each of the four implementations. An indicates that the solver failed to solve the 
problem instance. Such a failure occurred because the storage requirements for the 
explicit transformation kernel representation exceeded the memory available and the 
solver terminated gracefully. Solution times represent only the CPU time required 
to solve the problem, and exclude the initial problem input, the time required by 
the factorization heuristic to identify the factored row structure and the final output 
to record the solution. All figures are in CPU seconds. 

The formulations of three of the test problems were strongly influenced by the 
design of the target solver: the TANKER model possesses a strong GUB structure 
since it contains a set of schedule selection constraints for each ship (i.e., from a set 
of candidate schedules, select at most one). The AMM04H model is a multicom- 
modity capacitated transshipment problem and is thus best suited to a pure network 
factorization. The PHOENIXIO model design was shaped by the generalized net- 
work factorization paradigm. The nature of the factored row structures shown in 
Table (9.2) supports this assertion. In the TANKER model, the factored pure net- 
work rows are exactly the same as the the factored GUB rows, and the heuristic 
constructs a generalized network factored row set by identifying one additional row 
to be paired with each GUB row. The AMM04H model may be viewed as a GUB 
factorization with a relatively modest GUB set consisting of the joint capacitation 
constraints, or as a pure network factorization with a relatively large pure network 
(PN) factored row set. In the PHOENIXIO model, the dominant structure is clearly 
the generalized network row structure. 
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As we expect, the factorization is most successful when the model is wed to the 
solver. Although we are surprised to find the performance of (PNET) competitive 
with (GUB) on the TANKER model, (GUB) clearly dominates the X and GNET 
solvers. Similarly, (GNET) dominates on the PHOENIXIO model and (PNET) on 
the AMM04H model. We would be disappointed if the results were otherwise. 

We see evidence in Table (9.3) to suggest that the approach of using heuristics 
to automatically identify factored structure has its pitfalls. In a number of problem 
instances, although our heuristics identify significantly larger factored sets as we 
progress from the base system to (GNET), we see little improvement in computation 
times (INVEST, 0DSM1, TAM8). In fact, we see in the TANKER model that 
the temptation to confound the modelers intended GUB structure by specifying a 
generalized network factorization leads to disastrous consequences, even though this 
tactic doubles the size of the factored row set. Interestingly enough, when we apply 
the (GNET) solver to the (GUB) factored row set, we are able to solve the problem 
in 16.5 CPU seconds. This suggests that the “quality” of a row factorization is not 
completely specified by size alone. 

We are encouraged by the observation that the transition from the basic system 
to (GUB) to (PNET) seldom degrades solution times, even when doing so yields 
little gain in the number of additional factored rows. This seems to contradict 
popular folklore, which suggests that computation times worsen as the transition 
from unadorned Simplex to (GUB) to (PNET) is made unless the transition is 
accompanied by a substantial increase in the size of the factored portion of the 
model. In fact, computational testing reported by others is frequently limited to 
models in which the number of explicit rows is in the range of one to twenty (see 
e.g., Chen and Saigal [1977], Glover, Karney, Klingman and Russell [1978], Glover 
and Klingman [19S1] ). Our results are all the more remarkable given the the lack of 
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guidance from the modeler for the “intended” row factorization. We note, however, 
that results are mixed for the transition from (PNET) to (GNET), and it seems 
clear that the applicability of the (GNET) factorization is not as general as that of 
(PNET). 

Finally, we observe that in several models it is factorization that separates 
success from failure in solving the problem wdth a given allocation of computer 
resources. This fact alone may be reason to consider this approach in practice. 

Our second interest with respect to computation is in the memory requirements 
of our algorithms. Our design strategy allocates memory to the data structures 
which represent F n so that we may successfully represent the largest dimension of 
factored kernel that may possibly arise. Limited memory remains to store the repre- 
sentation of A^ 1 , the explicit transformation kernel. During the solution process, if 
we encounter a representation of A ^ which requires more memory than is available, 
failure occurs. We measure the size of the computer representation of A fj 1 and the 
amount of available computer storage in terms of the elements of Afj 1 that can be 
stored. The number of bytes per element varies according to the size of the problem 
(this has to do with the FORTRAN data types INTEGER*2 and INTEGER*4), but 
is generally 28 bytes per element. Table (9.4) lists for each problem instance/solver 
pair the maximum size of the A 7 / representation encountered during the solution, 
measured in units of number of elements. An asterisk (*) indicates that the number 
shown equals the maximum number of units of storage that were available, and thus 
failure occurred. 

We see that generally the representation of the maximum size of the explicit 
transformation kernel decreases as the generality of the factorization increases. Re- 
calling the definition of the explicit transformation kernel: 
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AJ = (E n -E n F l - l l Fn)-' 

this trend is as we would expect. As the generality of the factorization increases, 
we expect the size of the factored component to increase and the size of the ex- 
plicit component to decrease. Each potentially binding explicit row which may be 
converted to a factored row by adopting a more general factorization reduces the 
dimension of the resulting representation of An- Also, the density of the term 
—EiiF u 1 F\2 generally increases as the complexity of the structure of Fjj 1 increases. 
Assuming the dimension of Fn is k by k , the number of nonzeros in E^ 1 in the GUB 

factorization is k. In the pure network factorization, the number of nonzeros in Fn' 
2 

may be as large as y, and in the generalized network factorization, the number of 
nonzeros in E^ 1 may be as large as n 2 . We note that (GNET) again provides several 
exceptions to the general trend in Table (9.4). 

It is the dynamic nature of our factorization algorithms which marks this 
work as a departure from previous research. Table 9.5 illustrates the significance 
of this point. The first column lists the number of constraints which are binding 
at optimality, and the second column expresses this as a percentage of the total 
number of constraints in the problem instance. The results shown here are typical 
of real-world large-scale models. It is usually the case that many constraints are 
not binding at optimality, and there are computational advantages to be gained by 
exploiting this fact. 

Columns 3, 4 and 5 of Table 9.5 list for (GUB), (PNET) and (GNET) respec- 
tively the number of explicit constraints that are binding at optimality. Since in 
each implementation, binding factored constraints are handled more efficiently than 
binding explicit constraints, we see that the computational success of our dynamic 



155 



factorization algorithms is due to the fact that even in large model instances, we 
are able to limit our attention to a relatively small number of explicit constraints, 
usually on the order of a few hundred or less. While this is well beyond the size of 
previously reported implementations, our results show that it is quite manageable. 

It is useful to establish a criteria for comparing and contrasting the perfor- 
mance of the algorithms which accounts for both execution time and storage re- 
quirements. Although more sophisticated models could undoubtedly be developed, 
we offer a simple model which we feel captures the essential features we wish to 
consider. Define 



f(s,t) = s-t 

where s is the total computer storage (measured in megabytes) required to solve 
a problem instance (including program code, original problem data, tableau rep- 
resentation, factored kernel representation and explicit transformation kernel rep- 
resentation) and t is the execution time (measured in CPU seconds). f(s,t) is 
monotonically increasing in s and t, and in a crude fashion it captures the essential 
features of the way in which computer resources are often marketed commercially. 
We will use /(s, t ) as a measure of performance. Table (9.6) displays f(s,t) for each 
of the problem instance/solver combinations. indicates that the solver failed to 
solve the particular instance. 

We observe that the trend as the transitions are made from base system to 
(PNET) is a decline in f(s,t). It is apparent that (PNET) is a versatile imple- 
mentation, performing extremely well on models with highly favorable structure 
(AMM04H, KG4, GASPNA, GASPNC) and comparing favorably with (GUB) and 
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(GNET) on every other model. The performance of (GNET) is generally quite com- 
parable to that of (PNET). It appears to be about 15-20% slower than (PNET) 
when it is used to solve a pure network factorization (GASPNC, AMM04H). Our 
row elimination heuristics are usually effective in identifying a favorable factored 
structure, but the computational results on the TANKER model clearly illustrate 
the limitations of this approach. (GNET) apparently requires a more sophisticated 
and careful user than do (GUB) or (PNET), and in this sense it is perhaps a more 
specialized algorithm. The evidence shown here indicates that (PNET) is a strong 
candidate for use as a general implementation, and need not be viewed as a highly 
specialized implementation suitable only for rare model instances. For the models 
studied here, we have progressed well beyond the stage of solving instances with 
only a handful of explicit rows. 
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TABLE 9.2: Summary of Problem Suite Dimensions 





m 


n 


GTE 


960 


6,624 


INVEST 


1,338 


11,989 


TANKER 


83 


7,598 


GAS PN A 


6,848 


27,884 


GAS PN C 


3,794 


15,362 


GAS PN E 


1,184 


5,102 


GAS GN A 


6,848 


27,884 


GAS GN C 


3,974 


15,362 


GK 2 


3,819 


17,844 


GK 3 


5,728 


27,493 


GK 4 


7,636 


37,139 


ODSM1 


3,023 


11,568 


ODSM2 


594 


22,211 


TAM1 


91 


389 


TAM 2 


180 


1,204 


TAM 3 


269 


2,883 


TAM 4 


211 


1,327 


TAM5 


438 


10,969 


TAMS 


420 


6,104 


TAM12 


629 


17,793 


PHOENIXIO 


1,618 


6,884 


PHOENIX30 


4,305 


4,297 


AMMO 4H 


13,963 


83,497 



Ppn 


P GN 


NZ 


909 


917 


57,563 


1,101 


1,168 


36,829 


33 


66 


30,890 


5,934 


5,976 


36,702 


3,418 


3,420 


19,701 


877 


883 


7,355 


5,142 


5,976 


36,702 


3,084 


3,420 


19,701 


2,578 


2,585 


34,809 


3,867 


3,876 


54,289 


5,156 


5,167 


73,766 


540 


558 


21,532 


490 


490 


42,827 


28 


34 


1,212 


54 


66 


6,869 


81 


99 


21,356 


77 


98 


6,954 


132 


162 


93,964 


154 


196 


49,376 


231 


295 


164,947 


220 


1,153 


13,818 


303 


3,605 


16,441 


12,892 


12,892 


166,713 



Pgub 

909 

941 

33 

4,345 

2,658 

434 

4,484 

2,664 

1,265 

2,295 

2,428 

523 

490 

22 

42 

63 

59 

102 

118 

177 

206 

293 

1,071 
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TABLE 9.3: 


Solution 


Times 


in CPU 


Seconds 




X 


GUB 


PNET 


GNET 


GTE 


48.4 


33.0 


37.2 


38.3 


INVEST 


21.6 


17.1 


17.1 


17.0 


TANKER 


141.8 


20.5 


17.7 


43.8 


GAS PN A 


* 


* 


413.2 


527.7 


GAS PN C 


56.9 


53.7 


35.8 


31.3 


GAS PN E 


* 


146.8 


17.4 


20.7 


GAS GN A 


* 


* 


* 


521.6 


GAS GN C 


52.9 


50.9 


37.6 


33.7 


GK 2 


27.1 


22.8 


21.3 


19.4 


GI< 3 


62.4 


58.4 


50.4 


45.1 


GK 4 


* 


380.9 


182.2 


186.6 


ODSM1 


8.7 


12.1 


9.3 


8.8 


ODSM2 


34.1 


32.5 


32.8 


30.4 


TAM1 


0.7 


0.7 


0.7 


0.7 


TAM2 


2.5 


1.7 


2.2 


2.1 


TAM3 


17.8 


14.1 


14.3 


13.5 


TAM4 


1.2 


1.2 


1.5 


1.7 


TAM5 


317.1 


300.9 


307.8 


269.6 


TAM8 


93.3 


81.8 


83.8 


77.4 


TAM12 


660.9 


660.6 


624.1 


480.9 


PHOENIXIO 


68.6 


43 7 


30.4 


10.7 


PHOENIX30 


* 


* 


* 


495.1 


AMMO 4H 


NA f 


1,241.3 


253.9 


288.4 



• NA indicates problem instance not run. 

f This problem was run on an IBM 3081 K under the MVS operating system using 
VS FORTRAN 1.4.1 in a 32 megabyte virtual machine using an advanced starting 
solution. The solution time shown is adjusted to account for the approximate dif- 
ference in computing speed between the IBM 3081K and the IBM 3033/AP. 
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TABLE 9.4: Number of Elements in Explicit Transformation Kernel Rep- 
resentation of Optimality or at Failure 



X 



GTE 


11,571 


INVEST 


8,380 


TANKER 


2,256 


GAS PN A 


*132,995 


GAS PN C 


21,912 


GAS PN E 


*174,482 


GAS GN A 


*132,995 


GAS GN C 


17,966 


GK 2 


8,280 


GK 3 


13,991 


GK 4 


*110,334 


ODSM1 


1,322 


ODSM2 


418 


TAM1 


1,713 


TAM2 


3,042 


TAM 3 


10,430 


TAM4 


1,489 


TAM 5 


29,661 


TAM8 


19,545 


TAM12 


37,716 


PHOENIX 10 


73,737 


PHOENIX30 


*153,781 


AMMO 4H 


NA 



GUB 


PNET 


GNET 


221 


214 


196 


682 


654 


521 


415 


312 


81 


130,883 


35,395 


65,558 


14,765 


260 


364 


111,821 


4,808 


4,312 


130,883 


*127,780 


65,081 


8,017 


1,358 


608 


2,719 


221 


184 


7,771 


388 


310 


79,646 


6,532 


9,155 


861 


259 


421 


3 


3 


3 


1,231 


841 


1,165 


2,078 


1,444 


1,610 


8,326 


5,305 


6,525 


1,302 


867 


1,332 


23,790 


15,689 


18,737 


14,886 


8,351 


11,182 


28,935 


16,097 


20,486 


38,705 


21,879 


1,458 


151,746 


*148,348 


11,523 


122,085 


23 


23 
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TABLE 9.5: Summary of the Number of Binding Explicit Constraints at 
Optimality 





Total 

Binding 


Percent 


GUB 


PNET 


GNET 


GTE 


552 


57.5 


19 


18 


18 


INVEST 


758 


56.7 


199 


194 


162 


TANKER 


50 


60.2 


31 


30 


9 


GAS PN A 


3,051 


44.6 


* 


292 


349 


GAS PN C 


2,337 


61.6 


849 


66 


91 


GAS PN E 


897 


75.8 


572 


88 


81 


GAS GN A 


3,045 


44.5 


* 


★ 


352 


GAS GN C 


2,331 


58.7 


863 


402 


88 


GK 2 


1 ,950 


51.1 


1,011 


93 


90 


GI< 3 


2,912 


50.8 


1,360 


140 


141 


GK 4 


4,119 


53.9 


2,477 


363 


356 


0DSM1 


297 


9.8 


53 


49 


47 


0DSM2 


448 


75.4 


219 


215 


0 


TAM1 


46 


50.5 


44 


32 


35 


TAM2 


87 


48.3 


77 


51 


61 


TAM3 


148 


55.0 


131 


97 


107 


TAM4 


121 


57.3 


110 


59 


81 


TAM5 


252 


57.5 


228 


170 


186 


TAM 8 


276 


65.7 


238 


140 


174 


TAM 12 


413 


65.7 


355 


208 


252 


PHOENIXlO 


1,085 


67.1 


1,083 


1,082 


77 


PHOENIX30 


3,477 


80.8 


* 


* 


109 


AMMO 4H 


2,889 


20.7 


2,882 


7 


7 
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TABLE 9.6: f(s,t ) in Megabyte • seconds 



X GUB PNET GNET 



GTE 


84 


47 


54 


64 


INVEST 


25 


16 


19 


20 


TANKER 


112 


15 


13 


42 


GAS PN A 


* 


* 


926 


1,771 


GAS PN C 


73 


58 


27 


31 


GAS PN E 


* 


509 


9 


14 


GAS GN A 


* 


* 


* 


1,774 


GAS GN C 


62 


45 


30 


34 


GK 2 


27 


19 


18 


21 


GK 3 


90 


74 


60 


66 


GK 4 


* 


1,368 


309 


387 


ODSM1 


4 


6 


5 


7 


ODSM2 


49 


47 


48 


53 


TAM1 


0.1 


0.1 


0.1 


0.1 


TAM2 


1 


1 


1 


1 


TAM3 


14 


10 


9 


12 


TAM4 


0.1 


0.1 


0.1 


0.1 


TAM 5 


750 


662 


614 


622 


TAMS 


133 


106 


95 


110 


TAM12 


2,376 


2,213 


1,882 


1,632 


PHOENIXIO 


169 


65 


32 


7 


PHOENIX30 


* 


* 


* 


693 


AMMO 4H 


NA 


9,093 


933 


1194 
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X. CONCLUSIONS 



We have presented three dynamic row factorization algorithms for solving 
large-scale linear programs. Although each may be used to solve any LP instance, 
each is designed to exploit a particular model row structure: generalized upper 
bounds (GUB), pure network rows (PNET) or generalized network rows (GNET). 

Previous research by others generally suggests that specialized algorithms such 
as those presented here are useful only when the factored structure completely dom- 
inates the structure of the model instance. There are reports of algorithms for 
solving problems having a single unfactored (explicit) constraint (Hultz and Kling- 
man [1978], Klingman and Russell [1978]). When implementations are reported, 
problem suites are limited to instances having a very small number of explicit con- 
straints, typically in the range from one to twenty (Chen and Saigal [1977], Glover, 
Karney, Klingman and Russell [1978], Glover and Klingman [1981]). The consensus 
seems to be that such algorithms are appropriately viewed as specialized algorithms, 
useful only for solving very special problem instances. 

Our experience strongly refutes this view. We find the performances of our im- 
plementations of the dynamic factorization algorithms are competitive with that of 
a commercial-quality optimization system on every model instance we have tested. 
This is particularly remarkable for two reasons. First, our test suite consists of 
models developed by skilled modelers specifically to exploit the capabilities and 
characteristics of the solver with which our implementations are competing. Sec- 
ond, we must select the row factorizations without the benefit of guidance from the 
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modeler, relying instead on useful but imperfect heuristics. Despite these computa- 
tional handicaps, our tests show our implementations to be at least as efficient as a 
well-respected commercial-quality optimization system. 

Our development has stressed the similarity between the algorithms and the 
natural extension which leads from one to the next. This is in contrast to the devel- 
opment that has been reported for similar, non-dynamic algorithms (e.g., Dantzig 
and Van Slyke [1967], Klingman and Russell [1978] and Hultz and Klingman [1978]) 
in which the specifics of the individual algorithm obscure the generality of the ap- 
proach. The conceptual difference between our algorithms is seen to be largely 
isolated to the structure of a single algebraic entity, the factored kernel. By abstract- 
ing the structure of the factored kernel and concentrating on the general algorithm 
design, we demonstrate the versatility and flexibility of this approach. 

We are gratified to find that the modularity suggested by the algorithmic 
development can be realized in an implementation design. We succeed in developing 
a software suite wffiich displays a “single-system image”. The modularity of the 
algorithm allows the definition of an “abstract data type” (see, e.g., Aho, Hopcroft 
and Ullman [1974]) which isolates the data structures and update procedures for the 
factored kernel from the rest of the implementation. Each factorization is seamlessly 
integrated within the system design, presenting a single design image. 

The early 1980’s produced a great deal of research in the area of automatic 
identification of special structure in LP models (see, e.g., Gunawardane and Schrage 
[1977], Glover [1980], Schrage [1981], Brown, McBride and Wood [1985] and Bixby 
and Fourer [19S6]). We have incorporated the most useful of these ideas into our 
implementation, and we have what we believe to be the first implementation which 
supports the automatic identification of factored row sets. This capability may 
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be used to identify new factored structure or to validate or augment a modeler- 
provided recommendation. Our computational experience indicates that while this 
approach is not as promising as perhaps first envisioned, it is nonetheless a valuable 
tool. When faced with the choice of either solving an unfactored model instance 
or automatically identifying a factored structure and then using the corresponding 
solver, our results show that the latter is nearly always to be preferred. Our results 
seem to suggest, however, that in addition to quantity of factored rows, the issue 
of quality of factored rows exerts influence on the performance of the factorization 
algorithms. While not well understood, it is clear that the myopic approach of 
our heuristics is no substitute for the modeler’s guidance in identifying factored 
structure. 

Several areas suggest themselves for further research. Certainly additional fac- 
tored structures can be examined. For example, one approach for treating factored 
column structures (“complicating columns”) is to allow the partitioning of rows into 
the categories “factored” and “explicit” to vary as the algorithm progresses. That 
is, allow factored row set membership to be determined with respect to the column 
structure of the currently nonbasic variables rather than with respect to the column 
structure of all problem variables. While conceptually simple, such a generalization 
seems to present significant algorithmic challenges. 

General algorithms are sometimes useful in specialized contexts. For example, 
processing networks (Koene [1982]) are network models which allow proportional 
flow restrictions on the arcs entering or leaving some nodes. One formulation of such 
a model results in a pure or generalized network structure with a set of complicating 
columns. Chen and Engquist [1986] propose a primal partitioning algorithm for 
solving processing network problems. An alternate formulation yields a pure or 
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generalized structure with complicating rows, and we note that this is precisely the 
structure we seek for our network factorizations. 

The multicommodity capacitated transshipment problem (MCTP) has been 
the subject of much research over the years, and a number of specialized algorithms 
(see, e.g., Assad [1978] or Kennington [1978]) have been proposed to solve it. Adopt- 
ing a general perspective, MCTP may be viewed either as a GUB model or as a pure 
network model with side constraints, and either view might be preferred depending 
upon the dimensions of the particular instance under consideration. Our compu- 
tational experience indicates that the pure network factorization algorithm offers a 
powerful technique for solving MCTP. As an experiment, we customized our (PN) 
implementation to exploit the special structure of the side constraints in MCTP. It 
is interesting to note that in our scientific computing environment, we observed no 
difference in solution times between the customized version of (PN) and the original 
implementation. 

Finally, all the approaches we have considered assume the prior existence of a 
specific structure in the factored rows which in turn determines the structure of the 
factored kernel. An extension of this general approach is to relax the requirement for 
strict conformance to a specific structure. Instead, we might allow the factored row 
structure to be “nearly” homogeneous. For example, we may allow a small number 
of complicating columns to disturb what is otherwise a factored pure network row 
structure. We then expect the structure of the factored kernel to be dominated by 
that induced by the predominant row structure, with only occasional complications 
due to the exceptional row structure. We allow for this exceptional structure in 
the factored kernel by identifying it “on-the-fly” as the algorithm progresses, and 
treating it in an appropriate manner. This approach may be thought of as a hy- 
brid between the factored method developed here and dynamic basis triangulation 
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methods (see, e.g., Hellerman and Rarick [1971] and [1972], Saunders [1976] and 
McBride [1980]). 

Dynamic extrinsic factorization is subsumed if we activate functions in the 
update analogous to the secondary exchanges now employed. Essentially all that 
has to be done is ensure that successive factored components retain their stipulated 
special structure. In our estimation, this will only be justified in cases where the 
model structure is amenable, and quite likely will require some model-specific fea- 
tures to perform well on difficult models. We have limited our experimentation to 
those static extrinsic cases which are believed to be most useful. 
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